The Pyretechnics Fire Behavior Library
Copyright © 2023-2025 Spatial Informatics Group, LLC

import cython
import cython as cy
from functools import reduce
import numpy as np
if cython.compiled:
    from cython.cimports.numpy import ndarray
    from cython.cimports.pyretechnics.cy_types import pyidx
else:
    from numpy import ndarray
    from pyretechnics.py_types import pyidx

@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def is_pos_int(x: object) -> cy.bint:
    return isinstance(x, int) and x > 0


@cy.cfunc
def divide_evenly(dividend: cy.int, divisor: cy.int) -> cy.int:
    if divisor == 0:
        raise ValueError(f"{divisor} may not be zero.")
    else:
        quotient : cy.int = dividend // divisor
        remainder: cy.int = dividend % divisor
        if remainder == 0:
            return quotient
        else:
            raise ValueError(f"{dividend} must be an exact multiple of {divisor}.")


@cy.ccall
@cy.exceptval(check=False)
def to_positive_index_range(index_range: tuple[pyidx, pyidx]|None, axis_length: pyidx) -> tuple[pyidx, pyidx]:
    """
    Translate None and negative indices to positive indices.
    """
    if index_range is None:
        return (0, axis_length)
    else:
        start: pyidx|None = index_range[0]
        stop : pyidx|None = index_range[1]
        return (
            0 if start is None else axis_length + start if start < 0 else start,
            axis_length if stop is None else axis_length + stop if stop < 0 else stop
        )


@cy.cfunc
def stretch_array(old_array: ndarray, new_length: cy.int, repetitions: cy.float) -> ndarray:
    new_array: ndarray = np.zeros(new_length, dtype=old_array.dtype)
    i        : pyidx
    for i in range(new_length):
        new_array[i] = old_array[int(i / repetitions)]
    return new_array


@cy.cfunc
def maybe_repeat_array(maybe_array: ndarray, axis_repetitions: tuple[pyidx, cy.float]) -> ndarray:
    """
    Return a new array that is created by repeating the elements from the input
    array repetitions times along the specified array axis. Avoid allocating
    new memory if repetitions == 1 or if the repeated array axis has length 1.
    """
    axis       : pyidx    = axis_repetitions[0]
    repetitions: cy.float = axis_repetitions[1]
    old_array  : ndarray  = np.asarray(maybe_array)
    array_dims : pyidx    = old_array.ndim
    array_shape: list     = list(np.shape(old_array))
    axis_length: cy.int   = array_shape[axis]
    if repetitions == 1.0:
        # no repetitions necessary
        return old_array
    elif axis_length == 1:
        # broadcast single-element axis repetitions times
        array_shape[axis] = int(repetitions)
        return np.broadcast_to(old_array, array_shape)
    elif repetitions % 1.0 == 0.0:
        # repeat each element on the chosen axis repetitions times
        return np.repeat(old_array, int(repetitions), axis)
    elif axis == 0 and array_dims == 1:
        # populate a new 1D array of the expected length by translating its indices into the original array
        new_axis_length: cy.int = int(axis_length * repetitions)
        return stretch_array(old_array, new_axis_length, repetitions)
    elif axis == 0 and array_dims == 2:
        # populate a new 2D array of the expected length by translating its indices into the original array
        new_rows: cy.int = int(axis_length * repetitions)
        return np.stack([old_array[int(i / repetitions)] for i in range(new_rows)])
    elif axis == 1 and array_dims == 2:
        # populate a new 2D array of the expected length by translating its indices into the original array
        old_rows: cy.int = array_shape[0]
        new_cols: cy.int = int(axis_length * repetitions)
        return np.stack([stretch_array(old_array[i], new_cols, repetitions) for i in range(old_rows)])
    elif axis == 0 and array_dims == 3:
        # populate a new 3D array of the expected length by translating its indices into the original array
        new_bands: cy.int = int(axis_length * repetitions)
        return np.stack([old_array[int(b / repetitions)] for b in range(new_bands)])
    elif axis == 1 and array_dims == 3:
        # populate a new 2D array of the expected length by translating its indices into the original array
        old_bands: cy.int = array_shape[0]
        new_rows : cy.int = int(axis_length * repetitions)
        return np.stack([
            np.stack([
                old_array[b, int(i / repetitions)]
                for i in range(new_rows)
            ])
            for b in range(old_bands)
        ])
    elif axis == 2 and array_dims == 3:
        # populate a new 3D array of the expected length by translating its indices into the original array
        old_bands: cy.int = array_shape[0]
        old_rows : cy.int = array_shape[1]
        new_cols : cy.int = int(axis_length * repetitions)
        return np.stack([
            np.stack([
                stretch_array(old_array[b,i], new_cols, repetitions)
                for i in range(old_rows)
            ])
            for b in range(old_bands)
        ])
    else:
        raise ValueError("Floating point repetitions are only supported for 1D, 2D, and 3D arrays.")

@cy.cclass
class ISpaceTimeCube:
    @cy.cfunc
    @cy.exceptval(check=False)
    def get(self, t: pyidx, y: pyidx, x: pyidx) -> cy.float:
        pass

@cy.cclass
class SpaceTimeCube(ISpaceTimeCube):
    """
    Create an object that represents a 3D array with dimensions (T,Y,X) given by cube_shape.
    Internally, data is stored as a 3D Numpy array at the resolution of the provided base data.
    Whenever a point value or contiguous space-time region of values is requested, translate
    the given cube_shape coordinates into base coordinates, look up the values from the base data,
    expand them (if necessary) back into the cube_shape resolution, and return the resulting scalar
    value or array to the caller.
    """
    ndim         : cy.int
    size         : cy.ulonglong
    shape        : tuple[cy.int, cy.int, cy.int]
    base         : object
    t_repetitions: cy.float
    y_repetitions: cy.float
    x_repetitions: cy.float
    data         : cy.float[:,:,::1] # FIXME: Restore polymorphism for the underlying Numpy arrays


    def __init__(self, cube_shape: tuple[int, int, int], base: object) -> cy.void:
        """
        NOTE: The resolutions in cube_shape must be exact multiples of any existing dimensions
              in the base data. If base is not a Numpy float32 array, a new array will be allocated.
        """
        # Ensure that cube_shape contains 3 values or throw an error
        if len(cube_shape) != 3:
            raise ValueError("The cube_shape must contain exactly three values.")

        # Unpack the cube_shape values without type-checking
        cube_bands_: object = cube_shape[0]
        cube_rows_ : object = cube_shape[1]
        cube_cols_ : object = cube_shape[2]

        # Ensure that cube_shape only contains positive integers or throw an error
        if not(is_pos_int(cube_bands_) and is_pos_int(cube_rows_) and is_pos_int(cube_cols_)):
            raise ValueError("The cube_shape must contain only positive integers.")

        # Cast the cube_shape values as primitive ints
        cube_bands: cy.int = cube_bands_
        cube_rows : cy.int = cube_rows_
        cube_cols : cy.int = cube_cols_

        # Store the cube metadata for later
        self.ndim  = 3
        self.size  = cube_bands * cube_rows * cube_cols
        self.shape = (cube_bands, cube_rows, cube_cols)
        self.base  = base

        # Store the base data as a 3D array along with its axis repetitions
        base_dimensions: cy.int = np.ndim(base)

        if base_dimensions == 0:
            # 0D: Constant Input
            self.t_repetitions = cube_bands
            self.y_repetitions = cube_rows
            self.x_repetitions = cube_cols
            self.data          = np.asarray([[[base]]], dtype=np.float32)

        elif base_dimensions == 1:
            # 1D: Time-Series Input
            base_bands: cy.int = len(base)
            self.t_repetitions = float(cube_bands) / float(base_bands)
            self.y_repetitions = cube_rows
            self.x_repetitions = cube_cols
            # Ensure that the cube_shape is not smaller than the base shape
            if cube_bands < base_bands:
                raise ValueError("The cube_shape may not be smaller than the base shape.")
            # Warn if any repetitions are not whole numbers
            if self.t_repetitions % 1.0 != 0.0:
                print("WARNING: Input data's shape does not evenly divide the cube_shape."
                      + " Index lookups beyond the cube's edge may mistakenly return values without errors.")
            # Warn if base is not a Numpy float32 array
            if not(isinstance(base, np.ndarray)) or (base.dtype != np.float32):
                print("WARNING: Input data is not a Numpy float32 array. Data will be copied into SpaceTimeCube.",
                      flush=True)
            # Expand (base_bands) -> (base_bands,1,1)
            self.data = np.expand_dims(np.asarray(base, dtype=np.float32), axis=(1,2))

        elif base_dimensions == 2:
            # 2D: Spatial Input
            base_shape: tuple  = np.shape(base)
            base_rows : cy.int = base_shape[0]
            base_cols : cy.int = base_shape[1]
            self.t_repetitions = cube_bands
            self.y_repetitions = float(cube_rows) / float(base_rows)
            self.x_repetitions = float(cube_cols) / float(base_cols)
            # Ensure that the cube_shape is not smaller than the base shape
            if cube_rows < base_rows or cube_cols < base_cols:
                raise ValueError("The cube_shape may not be smaller than the base shape.")
            # Warn if any repetitions are not whole numbers
            if self.y_repetitions % 1.0 != 0.0 or self.x_repetitions % 1.0 != 0.0:
                print("WARNING: Input data's shape does not evenly divide the cube_shape."
                      + " Index lookups beyond the cube's edge may mistakenly return values without errors.")
            # Warn if base is not a Numpy float32 array
            if not(isinstance(base, np.ndarray)) or (base.dtype != np.float32):
                print("WARNING: Input data is not a Numpy float32 array. Data will be copied into SpaceTimeCube.",
                      flush=True)
            # Expand (base_rows,base_cols) -> (1,base_rows,base_cols)
            self.data = np.expand_dims(np.asarray(base, dtype=np.float32), axis=0)

        elif base_dimensions == 3:
            # 3D: Spatio-Temporal Input
            base_shape: tuple  = np.shape(base)
            base_bands: cy.int = base_shape[0]
            base_rows : cy.int = base_shape[1]
            base_cols : cy.int = base_shape[2]
            self.t_repetitions = float(cube_bands) / float(base_bands)
            self.y_repetitions = float(cube_rows) / float(base_rows)
            self.x_repetitions = float(cube_cols) / float(base_cols)
            # Ensure that the cube_shape is not smaller than the base shape
            if cube_bands < base_bands or cube_rows < base_rows or cube_cols < base_cols:
                raise ValueError("The cube_shape may not be smaller than the base shape.")
            # Warn if any repetitions are not whole numbers
            if self.t_repetitions % 1.0 != 0.0 or self.y_repetitions % 1.0 != 0.0 or self.x_repetitions % 1.0 != 0.0:
                print("WARNING: Input data's shape does not evenly divide the cube_shape."
                      + " Index lookups beyond the cube's edge may mistakenly return values without errors.")
            # Warn if base is not a Numpy float32 array
            if not(isinstance(base, np.ndarray)) or (base.dtype != np.float32):
                print("WARNING: Input data is not a Numpy float32 array. Data will be copied into SpaceTimeCube.",
                      flush=True)
            self.data = np.asarray(base, dtype=np.float32)

        else:
            # 4D+: Invalid Input
            raise ValueError("Invalid input: base must have 0-3 dimensions.")


    @cy.ccall
    @cy.exceptval(check=False)
    def get(self, t: pyidx, y: pyidx, x: pyidx) -> cy.float:
        """
        Return the scalar value at index (t,y,x) by translating these cube coordinates
        to base coordinates and looking up the value within the base data.

        NOTE: Indices may be negative.
        """
        # Select value by spatio-temporal coordinate
        base_t: pyidx = int(t / self.t_repetitions)
        base_y: pyidx = int(y / self.y_repetitions)
        base_x: pyidx = int(x / self.x_repetitions)
        return self.data[base_t, base_y, base_x]


    @cy.ccall
    def getTimeSeries(self, t_range: tuple[pyidx, pyidx]|None, y: pyidx, x: pyidx) -> ndarray:
        """
        Return the 1D array given by the slice (t_range,y,x) by translating these cube
        coordinates to base coordinates, looking up the array slice within the base data,
        and expanding it back to the cube_shape resolution.

        NOTE: Indices may be negative.
        NOTE: Range indices may include one or more None values and
              provide (inclusion, exclusion) semantics like Python array slice notation.
        """
        # Destructure the argument range
        t_range_updated : tuple[pyidx, pyidx] = to_positive_index_range(t_range, self.shape[0])
        t_start         : pyidx               = t_range_updated[0]
        t_stop_exclusive: pyidx               = t_range_updated[1]
        t_stop          : pyidx               = t_stop_exclusive - 1
        # Translate high-res coordinates to low-res coordinates
        t_start_chunk: pyidx = int(t_start / self.t_repetitions)
        t_stop_chunk : pyidx = int(t_stop  / self.t_repetitions)
        y_chunk      : pyidx = int(y       / self.y_repetitions)
        x_chunk      : pyidx = int(x       / self.x_repetitions)
        # Select the array slice that completely contains all low-res coordinates
        low_res_time: ndarray = np.asarray(self.data[t_start_chunk:(t_stop_chunk + 1),
                                                     y_chunk,
                                                     x_chunk])
        # Expand the low-res slice into a high-res slice
        high_res_time: ndarray = maybe_repeat_array(low_res_time, (0, self.t_repetitions))
        # Translate high-res global coordinates to high-res slice coordinates
        t_chunk_origin: pyidx = int(t_start_chunk * self.t_repetitions)
        t_start_idx   : pyidx = t_start - t_chunk_origin
        t_stop_idx    : pyidx = t_stop  - t_chunk_origin
        # Select the array slice that matches the high-res slice coordinates
        return high_res_time[t_start_idx:(t_stop_idx + 1)]


    @cy.ccall
    def getSpatialPlane(self,
                        t      : pyidx,
                        y_range: tuple[pyidx, pyidx]|None,
                        x_range: tuple[pyidx, pyidx]|None) -> ndarray:
        """
        Return the 2D array given by the slice (t,y_range,x_range) by translating these
        cube coordinates to base coordinates, looking up the array slice within the base
        data, and expanding it back to the cube_shape resolution.

        NOTE: Indices may be negative.
        NOTE: Range indices may include one or more None values and
              provide (inclusion, exclusion) semantics like Python array slice notation.
        """
        # Destructure the argument ranges
        y_range_updated : tuple[pyidx, pyidx] = to_positive_index_range(y_range, self.shape[1])
        x_range_updated : tuple[pyidx, pyidx] = to_positive_index_range(x_range, self.shape[2])
        y_start         : pyidx               = y_range_updated[0]
        y_stop_exclusive: pyidx               = y_range_updated[1]
        x_start         : pyidx               = x_range_updated[0]
        x_stop_exclusive: pyidx               = x_range_updated[1]
        y_stop          : pyidx               = y_stop_exclusive - 1
        x_stop          : pyidx               = x_stop_exclusive - 1
        # Translate high-res coordinates to low-res coordinates
        t_chunk      : pyidx = int(t       / self.t_repetitions)
        y_start_chunk: pyidx = int(y_start / self.y_repetitions)
        y_stop_chunk : pyidx = int(y_stop  / self.y_repetitions)
        x_start_chunk: pyidx = int(x_start / self.x_repetitions)
        x_stop_chunk : pyidx = int(x_stop  / self.x_repetitions)
        # Select the array slice that completely contains all low-res coordinates
        low_res_space: ndarray = np.asarray(self.data[t_chunk,
                                                      y_start_chunk:(y_stop_chunk + 1),
                                                      x_start_chunk:(x_stop_chunk + 1)])
        # Expand the low-res slice into a high-res slice
        high_res_space: ndarray = reduce(maybe_repeat_array,
                                         ((0, self.y_repetitions),
                                          (1, self.x_repetitions)),
                                         low_res_space)
        # Translate high-res global coordinates to high-res slice coordinates
        y_chunk_origin: pyidx = int(y_start_chunk * self.y_repetitions)
        x_chunk_origin: pyidx = int(x_start_chunk * self.x_repetitions)
        y_start_idx   : pyidx = y_start - y_chunk_origin
        y_stop_idx    : pyidx = y_stop  - y_chunk_origin
        x_start_idx   : pyidx = x_start - x_chunk_origin
        x_stop_idx    : pyidx = x_stop  - x_chunk_origin
        # Select the array slice that matches the high-res slice coordinates
        return high_res_space[y_start_idx:(y_stop_idx + 1),
                              x_start_idx:(x_stop_idx + 1)]


    @cy.ccall
    def getSubcube(self,
                   t_range: tuple[pyidx, pyidx]|None,
                   y_range: tuple[pyidx, pyidx]|None,
                   x_range: tuple[pyidx, pyidx]|None) -> ndarray:
        """
        Return the 3D array given by the slice (t_range,y_range,x_range) by translating
        these cube coordinates to base coordinates, looking up the array slice within the
        base data, and expanding it back to the cube_shape resolution.

        NOTE: Indices may be negative.
        NOTE: Range indices may include one or more None values and
              provide (inclusion, exclusion) semantics like Python array slice notation.
        """
        # Destructure the argument ranges
        t_range_updated : tuple[pyidx, pyidx] = to_positive_index_range(t_range, self.shape[0])
        y_range_updated : tuple[pyidx, pyidx] = to_positive_index_range(y_range, self.shape[1])
        x_range_updated : tuple[pyidx, pyidx] = to_positive_index_range(x_range, self.shape[2])
        t_start         : pyidx               = t_range_updated[0]
        t_stop_exclusive: pyidx               = t_range_updated[1]
        y_start         : pyidx               = y_range_updated[0]
        y_stop_exclusive: pyidx               = y_range_updated[1]
        x_start         : pyidx               = x_range_updated[0]
        x_stop_exclusive: pyidx               = x_range_updated[1]
        t_stop          : pyidx               = t_stop_exclusive - 1
        y_stop          : pyidx               = y_stop_exclusive - 1
        x_stop          : pyidx               = x_stop_exclusive - 1
        # Translate high-res coordinates to low-res coordinates
        t_start_chunk: pyidx = int(t_start / self.t_repetitions)
        t_stop_chunk : pyidx = int(t_stop  / self.t_repetitions)
        y_start_chunk: pyidx = int(y_start / self.y_repetitions)
        y_stop_chunk : pyidx = int(y_stop  / self.y_repetitions)
        x_start_chunk: pyidx = int(x_start / self.x_repetitions)
        x_stop_chunk : pyidx = int(x_stop  / self.x_repetitions)
        # Select the array slice that completely contains all low-res coordinates
        low_res_cube: ndarray = np.asarray(self.data[t_start_chunk:(t_stop_chunk + 1),
                                                     y_start_chunk:(y_stop_chunk + 1),
                                                     x_start_chunk:(x_stop_chunk + 1)])
        # Expand the low-res slice into a high-res slice
        high_res_cube: ndarray = reduce(maybe_repeat_array,
                                        ((0, self.t_repetitions),
                                         (1, self.y_repetitions),
                                         (2, self.x_repetitions)),
                                        low_res_cube)
        # Translate high-res global coordinates to high-res slice coordinates
        t_chunk_origin: pyidx = int(t_start_chunk * self.t_repetitions)
        y_chunk_origin: pyidx = int(y_start_chunk * self.y_repetitions)
        x_chunk_origin: pyidx = int(x_start_chunk * self.x_repetitions)
        t_start_idx   : pyidx = t_start - t_chunk_origin
        t_stop_idx    : pyidx = t_stop  - t_chunk_origin
        y_start_idx   : pyidx = y_start - y_chunk_origin
        y_stop_idx    : pyidx = y_stop  - y_chunk_origin
        x_start_idx   : pyidx = x_start - x_chunk_origin
        x_stop_idx    : pyidx = x_stop  - x_chunk_origin
        # Select the array slice that matches the high-res slice coordinates
        return high_res_cube[t_start_idx:(t_stop_idx + 1),
                             y_start_idx:(y_stop_idx + 1),
                             x_start_idx:(x_stop_idx + 1)]


    def __getFullyRealizedCube(self):
        """
        Return the 3D array created by expanding the base data to the cube_shape resolution.
        Wherever possible, Numpy broadcasting is used to avoid memory allocation along
        constant array dimensions.
        """
        base_dimensions = np.ndim(self.base)

        if base_dimensions == 0:
            # 0D: Constant Input
            # Broadcast (0,0,0) -> (t,y,x)
            return np.broadcast_to(self.data, self.shape)

        elif base_dimensions == 1:
            # 1D: Time-Series Input
            # Repeat (t0,1,1) -> (t,1,1)
            repeated_array = maybe_repeat_array(np.asarray(self.data), (0, self.t_repetitions))
            # Broadcast (t,1,1) -> (t,y,x)
            return np.broadcast_to(repeated_array, self.shape)

        elif base_dimensions == 2:
            # 2D: Spatial Input
            # Repeat (1,y0,x0) -> (1,y,x)
            repeated_array = reduce(maybe_repeat_array,
                                    ((1, self.y_repetitions),
                                     (2, self.x_repetitions)),
                                    np.asarray(self.data))
            # Broadcast (1,y,x) -> (t,y,x)
            return np.broadcast_to(repeated_array, self.shape)

        else:
            # 3D: Spatio-Temporal Input
            # Repeat (t0,y0,x0) -> (t,y,x)
            return reduce(maybe_repeat_array,
                          ((0, self.t_repetitions),
                           (1, self.y_repetitions),
                           (2, self.x_repetitions)),
                          np.asarray(self.data))


    def getFullyRealizedCube(self, cache=False):
        """
        Return the 3D array created by expanding the base data to the cube_shape resolution.
        Wherever possible, Numpy broadcasting is used to avoid memory allocation along
        constant array dimensions. When cache == True, this expanded 3D array is cached
        within the SpaceTimeCube object for future immediate retrieval.
        """
        if hasattr(self, "cube"):
            return self.cube
        else:
            cube = self.__getFullyRealizedCube()
            if cache is True:
                self.cube = cube
            return cube


    def releaseFullyRealizedCube(self):
        """
        Deletes the cached fully realized cube if it exists.
        """
        if hasattr(self, "cube"):
            delattr(self, "cube")

@cy.cclass
class LazySpaceTimeCube(ISpaceTimeCube):
    """
    Create an object that represents a 3D array with dimensions (T,Y,X) given by cube_shape.
    Internally, data is stored as an initially empty 3D array of SpaceTimeCube objects.
    Whenever a point value or contiguous space-time region of values is requested, identify
    which SpaceTimeCubes contain the requested coordinates, load them into the cache array
    by calling load_subcube for any that are not already present, request the values from
    these SpaceTimeCubes, combine them together if necessary, and return the resulting scalar
    value or array to the caller.
    """
    ndim         : cy.int
    size         : cy.ulonglong
    shape        : tuple[cy.int, cy.int, cy.int]
    subcube_shape: tuple[cy.int, cy.int, cy.int]
    cache_shape  : tuple[cy.int, cy.int, cy.int]
    cache        : ndarray
    load_subcube : object


    def __init__(self,
                 cube_shape   : tuple[int, int, int],
                 subcube_shape: tuple[int, int, int],
                 load_subcube : object) -> cy.void:
        """
        NOTE: The resolutions in cube_shape must be exact multiples of those in subcube_shape.
        """
        # Ensure that cube_shape and subcube_shape both contain 3 values or throw an error
        if len(cube_shape) != 3:
            raise ValueError("The cube_shape must contain exactly three values.")

        if len(subcube_shape) != 3:
            raise ValueError("The subcube_shape must contain exactly three values.")

        # Unpack the cube_shape values without type-checking
        cube_bands_: object = cube_shape[0]
        cube_rows_ : object = cube_shape[1]
        cube_cols_ : object = cube_shape[2]

        # Unpack the subcube_shape values without type-checking
        subcube_bands_: object = subcube_shape[0]
        subcube_rows_ : object = subcube_shape[1]
        subcube_cols_ : object = subcube_shape[2]

        # Ensure that cube_shape and subcube_shape only contain positive integers or throw an error
        if not(is_pos_int(cube_bands_) and is_pos_int(cube_rows_) and is_pos_int(cube_cols_)):
            raise ValueError("The cube_shape must contain only positive integers.")

        if not(is_pos_int(subcube_bands_) and is_pos_int(subcube_rows_) and is_pos_int(subcube_cols_)):
            raise ValueError("The subcube_shape must contain only positive integers.")

        # Cast the cube_shape values as primitive ints
        cube_bands: cy.int = cube_bands_
        cube_rows : cy.int = cube_rows_
        cube_cols : cy.int = cube_cols_

        # Cast the subcube_shape values as primitive ints
        subcube_bands: cy.int = subcube_bands_
        subcube_rows : cy.int = subcube_rows_
        subcube_cols : cy.int = subcube_cols_

        # Ensure that cube_shape is divided evenly by subcube_shape or throw an error
        cache_bands: cy.int = divide_evenly(cube_bands, subcube_bands)
        cache_rows : cy.int = divide_evenly(cube_rows, subcube_rows)
        cache_cols : cy.int = divide_evenly(cube_cols, subcube_cols)

        # Store the cube metadata, subcube_shape, cache_shape, cache, and load_subcube functions for later
        self.ndim          = 3
        self.size          = cube_bands * cube_rows * cube_cols
        self.shape         = (cube_bands, cube_rows, cube_cols)
        self.subcube_shape = (subcube_bands, subcube_rows, subcube_cols)
        self.cache_shape   = (cache_bands, cache_rows, cache_cols)
        self.cache         = np.empty(self.cache_shape, dtype=object)
        self.load_subcube  = load_subcube


    @cy.cfunc
    def __getOrLoadSubcube(self, cache_t: pyidx, cache_y: pyidx, cache_x: pyidx) -> SpaceTimeCube:
        """
        Return the SpaceTimeCube stored at self.cache[cache_t, cache_y, cache_x] if it
        has already been loaded. Otherwise, call self.load_subcube to load it, store
        it in self.cache, and return it.
        """
        subcube: SpaceTimeCube = cy.cast(SpaceTimeCube, self.cache[cache_t, cache_y, cache_x])
        if subcube:
            return subcube
        else:
            subcube = self.load_subcube((cache_t, cache_y, cache_x), self.subcube_shape)
            self.cache[cache_t, cache_y, cache_x] = subcube
            return subcube


    @cy.ccall
    @cy.exceptval(check=False)
    def get(self, t: pyidx, y: pyidx, x: pyidx) -> cy.float:
        """
        Return the scalar value at index (t,y,x) by translating these cube coordinates
        to cache and subcube coordinates, loading the matching subcube into the cache grid
        if not already present, and looking up the value within this subcube.

        NOTE: Indices may be negative provided that your load_subcube function can handle
              negative indices in its cache_index argument.
        """
        # Grab the subcube_shape tuple
        subcube_shape: tuple[cy.int, cy.int, cy.int] = self.subcube_shape

        # Unpack the subcube_shape values
        subcube_bands: cy.int = subcube_shape[0]
        subcube_rows : cy.int = subcube_shape[1]
        subcube_cols : cy.int = subcube_shape[2]

        # Calculate the cache index
        cache_t: pyidx = t // subcube_bands
        cache_y: pyidx = y // subcube_rows
        cache_x: pyidx = x // subcube_cols

        # Calculate the subcube index
        subcube_t: pyidx = t % subcube_bands
        subcube_y: pyidx = y % subcube_rows
        subcube_x: pyidx = x % subcube_cols

        # Fetch the subcube from the cache
        subcube: SpaceTimeCube = self.__getOrLoadSubcube(cache_t, cache_y, cache_x)

        # Look up the scalar value in the subcube at the subcube index
        return subcube.get(subcube_t, subcube_y, subcube_x)


    def getTimeSeries(self, t_range, y, x):
        """
        Return the 1D array given by the slice (t_range,y,x) by translating these cube
        coordinates to cache and subcube coordinates, loading the matching subcubes into
        the cache grid if not already present, looking up the array slices within each
        subcube, and merging them together into a single 1D array.

        NOTE: Indices may be negative provided that your load_subcube function can handle
              negative indices in its cache_index argument.
        NOTE: Range indices may include one or more None values and
              provide (inclusion, exclusion) semantics like Python array slice notation.
        """
        # Destructure the argument range
        (t_start, t_stop_exclusive) = to_positive_index_range(t_range, self.shape[0])
        t_stop = t_stop_exclusive - 1
        # Translate high-res coordinates to cache and subcube coordinates
        (subcube_bands, subcube_rows, subcube_cols) = self.subcube_shape
        (cache_t_start, subcube_t_start) = divmod(t_start, subcube_bands)
        (cache_t_stop,  subcube_t_stop)  = divmod(t_stop,  subcube_bands)
        (cache_y,       subcube_y)       = divmod(y,       subcube_rows)
        (cache_x,       subcube_x)       = divmod(x,       subcube_cols)
        # Load, expand, and combine subcubes
        return np.concatenate(
            [self.__getOrLoadSubcube(cache_t,
                                     cache_y,
                                     cache_x
                                    ).getTimeSeries(
                                        (subcube_t_start    if cache_t == cache_t_start else 0,
                                         subcube_t_stop + 1 if cache_t == cache_t_stop  else subcube_bands),
                                        subcube_y,
                                        subcube_x
                                    )
             for cache_t in range(cache_t_start, cache_t_stop + 1)]
        )


    def getSpatialPlane(self, t, y_range, x_range):
        """
        Return the 2D array given by the slice (t,y_range,x_range) by translating these
        cube coordinates to cache and subcube coordinates, loading the matching subcubes
        into the cache grid if not already present, looking up the array slices within each
        subcube, and merging them together into a single 2D array.

        NOTE: Indices may be negative provided that your load_subcube function can handle
              negative indices in its cache_index argument.
        NOTE: Range indices may include one or more None values and
              provide (inclusion, exclusion) semantics like Python array slice notation.
        """
        # Destructure the argument ranges
        (y_start, y_stop_exclusive) = to_positive_index_range(y_range, self.shape[1])
        (x_start, x_stop_exclusive) = to_positive_index_range(x_range, self.shape[2])
        y_stop = y_stop_exclusive - 1
        x_stop = x_stop_exclusive - 1
        # Translate high-res coordinates to cache and subcube coordinates
        (subcube_bands, subcube_rows, subcube_cols) = self.subcube_shape
        (cache_t,       subcube_t)       = divmod(t,       subcube_bands)
        (cache_y_start, subcube_y_start) = divmod(y_start, subcube_rows)
        (cache_y_stop,  subcube_y_stop)  = divmod(y_stop,  subcube_rows)
        (cache_x_start, subcube_x_start) = divmod(x_start, subcube_cols)
        (cache_x_stop,  subcube_x_stop)  = divmod(x_stop,  subcube_cols)
        # Load, expand, and combine subcubes
        return np.block(
            [[self.__getOrLoadSubcube(cache_t,
                                      cache_y,
                                      cache_x
                                      ).getSpatialPlane(
                                          subcube_t,
                                          (subcube_y_start    if cache_y == cache_y_start else 0,
                                           subcube_y_stop + 1 if cache_y == cache_y_stop  else subcube_rows),
                                          (subcube_x_start    if cache_x == cache_x_start else 0,
                                           subcube_x_stop + 1 if cache_x == cache_x_stop  else subcube_cols)
                                      )
              for cache_x in range(cache_x_start, cache_x_stop + 1)]
             for cache_y in range(cache_y_start, cache_y_stop + 1)]
        )


    def getSubcube(self, t_range, y_range, x_range):
        """
        Return the 3D array given by the slice (t_range,y_range,x_range) by translating
        these cube coordinates to cache and subcube coordinates, loading the matching
        subcubes into the cache grid if not already present, looking up the array slices
        within each subcube, and merging them together into a single 3D array.

        NOTE: Indices may be negative provided that your load_subcube function can handle
              negative indices in its cache_index argument.
        NOTE: Range indices may include one or more None values and
              provide (inclusion, exclusion) semantics like Python array slice notation.
        """
        # Destructure the argument ranges
        (t_start, t_stop_exclusive) = to_positive_index_range(t_range, self.shape[0])
        (y_start, y_stop_exclusive) = to_positive_index_range(y_range, self.shape[1])
        (x_start, x_stop_exclusive) = to_positive_index_range(x_range, self.shape[2])
        t_stop = t_stop_exclusive - 1
        y_stop = y_stop_exclusive - 1
        x_stop = x_stop_exclusive - 1
        # Translate high-res coordinates to cache and subcube coordinates
        (subcube_bands, subcube_rows, subcube_cols) = self.subcube_shape
        (cache_t_start, subcube_t_start) = divmod(t_start, subcube_bands)
        (cache_t_stop,  subcube_t_stop)  = divmod(t_stop,  subcube_bands)
        (cache_y_start, subcube_y_start) = divmod(y_start, subcube_rows)
        (cache_y_stop,  subcube_y_stop)  = divmod(y_stop,  subcube_rows)
        (cache_x_start, subcube_x_start) = divmod(x_start, subcube_cols)
        (cache_x_stop,  subcube_x_stop)  = divmod(x_stop,  subcube_cols)
        # Load, expand, and combine subcubes
        return np.block(
            [[[self.__getOrLoadSubcube(cache_t,
                                       cache_y,
                                       cache_x
                                       ).getSubcube(
                                           (subcube_t_start    if cache_t == cache_t_start else 0,
                                            subcube_t_stop + 1 if cache_t == cache_t_stop  else subcube_bands),
                                           (subcube_y_start    if cache_y == cache_y_start else 0,
                                            subcube_y_stop + 1 if cache_y == cache_y_stop  else subcube_rows),
                                           (subcube_x_start    if cache_x == cache_x_start else 0,
                                            subcube_x_stop + 1 if cache_x == cache_x_stop  else subcube_cols)
                                       )
               for cache_x in range(cache_x_start, cache_x_stop + 1)]
              for cache_y in range(cache_y_start, cache_y_stop + 1)]
             for cache_t in range(cache_t_start, cache_t_stop + 1)]
        )


    def getFullyRealizedCube(self, cache=False):
        raise ValueError("getFullyRealizedCube is not implemented for LazySpaceTimeCube.\n"
                         + "You probably don't want to do this anyway.")


    def releaseFullyRealizedCube(self):
        raise ValueError("releaseFullyRealizedCube is not implemented for LazySpaceTimeCube.\n"
                         + "You probably don't want to do this anyway.")

import cython
import cython as cy
if cython.compiled:
    from cython.cimports.libc.math import exp
    from cython.cimports.pyretechnics.cy_types import fcatarr, fclaarr, CompactFuelModel, FuelModel
else:
    from math import exp
    from pyretechnics.py_types import fcatarr, fclaarr, CompactFuelModel, FuelModel

# TODO: OPTIM Replace this dictionary with something more efficient
# Lookup table including entries for each of the Anderson 13 and Scott & Burgan 40 fuel models.
#
# The fields have the following meanings:
#   {
#     fuel_model_number: (delta,
#                         M_x_dead,
#                         h,
#                         w_o_dead_1hr,
#                         w_o_dead_10hr,
#                         w_o_dead_100hr,
#                         w_o_live_herbaceous,
#                         w_o_live_woody,
#                         sigma_dead_1hr,
#                         sigma_dead_10hr,
#                         sigma_dead_100hr,
#                         sigma_live_herbaceous,
#                         sigma_live_woody), # name
#   }
compact_fuel_model_table = cy.declare(dict, { # dict[int, CompactFuelModel]
    # Anderson 13:
    # Grass and Grass-dominated (short-grass,timber-grass-and-understory,tall-grass)
    1  : (1.0, 12.0, 8.0, 0.0340, 0.0000, 0.0000, 0.0000, 0.0000, 3500.0,   0.0,  0.0,    0.0,    0.0), # R01
    2  : (1.0, 15.0, 8.0, 0.0920, 0.0460, 0.0230, 0.0230, 0.0000, 3000.0, 109.0, 30.0, 1500.0,    0.0), # R02
    3  : (2.5, 25.0, 8.0, 0.1380, 0.0000, 0.0000, 0.0000, 0.0000, 1500.0,   0.0,  0.0,    0.0,    0.0), # R03
    # Chaparral and Shrubfields (chaparral,brush,dormant-brush-hardwood-slash,southern-rough)
    4  : (6.0, 20.0, 8.0, 0.2300, 0.1840, 0.0920, 0.2300, 0.0000, 2000.0, 109.0, 30.0, 1500.0,    0.0), # R04
    5  : (2.0, 20.0, 8.0, 0.0460, 0.0230, 0.0000, 0.0920, 0.0000, 2000.0, 109.0,  0.0, 1500.0,    0.0), # R05
    6  : (2.5, 25.0, 8.0, 0.0690, 0.1150, 0.0920, 0.0000, 0.0000, 1750.0, 109.0, 30.0,    0.0,    0.0), # R06
    7  : (2.5, 40.0, 8.0, 0.0520, 0.0860, 0.0690, 0.0170, 0.0000, 1750.0, 109.0, 30.0, 1550.0,    0.0), # R07
    # Timber Litter (closed-timber-litter,hardwood-litter,timber-litter-and-understory)
    8  : (0.2, 30.0, 8.0, 0.0690, 0.0460, 0.1150, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # R08
    9  : (0.2, 25.0, 8.0, 0.1340, 0.0190, 0.0070, 0.0000, 0.0000, 2500.0, 109.0, 30.0,    0.0,    0.0), # R09
    10 : (1.0, 25.0, 8.0, 0.1380, 0.0920, 0.2300, 0.0920, 0.0000, 2000.0, 109.0, 30.0, 1500.0,    0.0), # R10
    # Logging Slash (light-logging-slash,medium-logging-slash,heavy-logging-slash)
    11 : (1.0, 15.0, 8.0, 0.0690, 0.2070, 0.2530, 0.0000, 0.0000, 1500.0, 109.0, 30.0,    0.0,    0.0), # R11
    12 : (2.3, 20.0, 8.0, 0.1840, 0.6440, 0.7590, 0.0000, 0.0000, 1500.0, 109.0, 30.0,    0.0,    0.0), # R12
    13 : (3.0, 25.0, 8.0, 0.3220, 1.0580, 1.2880, 0.0000, 0.0000, 1500.0, 109.0, 30.0,    0.0,    0.0), # R13
    # Nonburnable (NB)
    91 : (0.0,  0.0, 0.0, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,    0.0,   0.0,  0.0,    0.0,    0.0), # NB1
    92 : (0.0,  0.0, 0.0, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,    0.0,   0.0,  0.0,    0.0,    0.0), # NB2
    93 : (0.0,  0.0, 0.0, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,    0.0,   0.0,  0.0,    0.0,    0.0), # NB3
    98 : (0.0,  0.0, 0.0, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,    0.0,   0.0,  0.0,    0.0,    0.0), # NB4
    99 : (0.0,  0.0, 0.0, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,    0.0,   0.0,  0.0,    0.0,    0.0), # NB5
    # Scott & Burgan 40:
    # Grass (GR)
    101: (0.4, 15.0, 8.0, 0.0046, 0.0000, 0.0000, 0.0138, 0.0000, 2200.0, 109.0, 30.0, 2000.0,    0.0), # GR1
    102: (1.0, 15.0, 8.0, 0.0046, 0.0000, 0.0000, 0.0459, 0.0000, 2000.0, 109.0, 30.0, 1800.0,    0.0), # GR2
    103: (2.0, 30.0, 8.0, 0.0046, 0.0184, 0.0000, 0.0689, 0.0000, 1500.0, 109.0, 30.0, 1300.0,    0.0), # GR3
    104: (2.0, 15.0, 8.0, 0.0115, 0.0000, 0.0000, 0.0872, 0.0000, 2000.0, 109.0, 30.0, 1800.0,    0.0), # GR4
    105: (1.5, 40.0, 8.0, 0.0184, 0.0000, 0.0000, 0.1148, 0.0000, 1800.0, 109.0, 30.0, 1600.0,    0.0), # GR5
    106: (1.5, 40.0, 9.0, 0.0046, 0.0000, 0.0000, 0.1561, 0.0000, 2200.0, 109.0, 30.0, 2000.0,    0.0), # GR6
    107: (3.0, 15.0, 8.0, 0.0459, 0.0000, 0.0000, 0.2479, 0.0000, 2000.0, 109.0, 30.0, 1800.0,    0.0), # GR7
    108: (4.0, 30.0, 8.0, 0.0230, 0.0459, 0.0000, 0.3352, 0.0000, 1500.0, 109.0, 30.0, 1300.0,    0.0), # GR8
    109: (5.0, 40.0, 8.0, 0.0459, 0.0459, 0.0000, 0.4132, 0.0000, 1800.0, 109.0, 30.0, 1600.0,    0.0), # GR9
    # Grass-Shrub (GS)
    121: (0.9, 15.0, 8.0, 0.0092, 0.0000, 0.0000, 0.0230, 0.0298, 2000.0, 109.0, 30.0, 1800.0, 1800.0), # GS1
    122: (1.5, 15.0, 8.0, 0.0230, 0.0230, 0.0000, 0.0275, 0.0459, 2000.0, 109.0, 30.0, 1800.0, 1800.0), # GS2
    123: (1.8, 40.0, 8.0, 0.0138, 0.0115, 0.0000, 0.0666, 0.0574, 1800.0, 109.0, 30.0, 1600.0, 1600.0), # GS3
    124: (2.1, 40.0, 8.0, 0.0872, 0.0138, 0.0046, 0.1561, 0.3260, 1800.0, 109.0, 30.0, 1600.0, 1600.0), # GS4
    # Shrub (SH)
    141: (1.0, 15.0, 8.0, 0.0115, 0.0115, 0.0000, 0.0069, 0.0597, 2000.0, 109.0, 30.0, 1800.0, 1600.0), # SH1
    142: (1.0, 15.0, 8.0, 0.0620, 0.1102, 0.0344, 0.0000, 0.1768, 2000.0, 109.0, 30.0,    0.0, 1600.0), # SH2
    143: (2.4, 40.0, 8.0, 0.0207, 0.1377, 0.0000, 0.0000, 0.2847, 1600.0, 109.0, 30.0,    0.0, 1400.0), # SH3
    144: (3.0, 30.0, 8.0, 0.0390, 0.0528, 0.0092, 0.0000, 0.1171, 2000.0, 109.0, 30.0, 1800.0, 1600.0), # SH4
    145: (6.0, 15.0, 8.0, 0.1653, 0.0964, 0.0000, 0.0000, 0.1331,  750.0, 109.0, 30.0,    0.0, 1600.0), # SH5
    146: (2.0, 30.0, 8.0, 0.1331, 0.0666, 0.0000, 0.0000, 0.0643,  750.0, 109.0, 30.0,    0.0, 1600.0), # SH6
    147: (6.0, 15.0, 8.0, 0.1607, 0.2433, 0.1010, 0.0000, 0.1561,  750.0, 109.0, 30.0,    0.0, 1600.0), # SH7
    148: (3.0, 40.0, 8.0, 0.0941, 0.1561, 0.0390, 0.0000, 0.1997,  750.0, 109.0, 30.0,    0.0, 1600.0), # SH8
    149: (4.4, 40.0, 8.0, 0.2066, 0.1125, 0.0000, 0.0712, 0.3214,  750.0, 109.0, 30.0, 1800.0, 1500.0), # SH9
    # Timber-Understory (TU)
    161: (0.6, 20.0, 8.0, 0.0092, 0.0413, 0.0689, 0.0092, 0.0413, 2000.0, 109.0, 30.0, 1800.0, 1600.0), # TU1
    162: (1.0, 30.0, 8.0, 0.0436, 0.0826, 0.0574, 0.0000, 0.0092, 2000.0, 109.0, 30.0,    0.0, 1600.0), # TU2
    163: (1.3, 30.0, 8.0, 0.0505, 0.0069, 0.0115, 0.0298, 0.0505, 1800.0, 109.0, 30.0, 1600.0, 1400.0), # TU3
    164: (0.5, 12.0, 8.0, 0.2066, 0.0000, 0.0000, 0.0000, 0.0918, 2300.0, 109.0, 30.0,    0.0, 2000.0), # TU4
    165: (1.0, 25.0, 8.0, 0.1837, 0.1837, 0.1377, 0.0000, 0.1377, 1500.0, 109.0, 30.0,    0.0,  750.0), # TU5
    # Timber Litter (TL)
    181: (0.2, 30.0, 8.0, 0.0459, 0.1010, 0.1653, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # TL1
    182: (0.2, 25.0, 8.0, 0.0643, 0.1056, 0.1010, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # TL2
    183: (0.3, 20.0, 8.0, 0.0230, 0.1010, 0.1286, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # TL3
    184: (0.4, 25.0, 8.0, 0.0230, 0.0689, 0.1928, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # TL4
    185: (0.6, 25.0, 8.0, 0.0528, 0.1148, 0.2020, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0, 1600.0), # TL5
    186: (0.3, 25.0, 8.0, 0.1102, 0.0551, 0.0551, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # TL6
    187: (0.4, 25.0, 8.0, 0.0138, 0.0643, 0.3719, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # TL7
    188: (0.3, 35.0, 8.0, 0.2663, 0.0643, 0.0505, 0.0000, 0.0000, 1800.0, 109.0, 30.0,    0.0,    0.0), # TL8
    189: (0.6, 35.0, 8.0, 0.3053, 0.1515, 0.1905, 0.0000, 0.0000, 1800.0, 109.0, 30.0,    0.0, 1600.0), # TL9
    # Slash-Blowdown (SB)
    201: (1.0, 25.0, 8.0, 0.0689, 0.1377, 0.5051, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # SB1
    202: (1.0, 25.0, 8.0, 0.2066, 0.1951, 0.1837, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # SB2
    203: (1.2, 25.0, 8.0, 0.2525, 0.1263, 0.1377, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # SB3
    204: (2.7, 25.0, 8.0, 0.2410, 0.1607, 0.2410, 0.0000, 0.0000, 2000.0, 109.0, 30.0,    0.0,    0.0), # SB4
})

# NOTE: We use this variable for X > 0.0 comparisons to account for floating point precision issues.
almost_zero = cy.declare(cy.float, 1e-6)


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def is_burnable_fuel_model_number(fuel_model_number: cy.int) -> cy.bint:
    return not (91 <= fuel_model_number <= 99)


@cy.cfunc
@cy.exceptval(check=False)
def compute_exp_A_sigma(A: cy.float, sigma_ij: cy.float) -> cy.float:
    if sigma_ij > almost_zero:
        return exp(A / sigma_ij)
    else:
        return 0.0


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def compute_firemod_size_class(sigma_i: cy.float) -> cy.float:
    return (
        1.0 if (sigma_i >= 1200.0)
        else 2.0 if (sigma_i >= 192.0)
        else 3.0 if (sigma_i >= 96.0)
        else 4.0 if (sigma_i >= 48.0)
        else 5.0 if (sigma_i >= 16.0)
        else 6.0
    )


@cy.cfunc
@cy.exceptval(check=False)
def expand_compact_fuel_model(fuel_model_number: cy.int) -> FuelModel:
    # Look up the CompactFuelModel by fuel_model_number
    cfm: CompactFuelModel = compact_fuel_model_table[fuel_model_number]
    # Unpack the CompactFuelModel values
    delta                : cy.float = cfm[0]
    M_x_dead             : cy.float = cfm[1]
    h                    : cy.float = cfm[2]
    w_o_dead_1hr         : cy.float = cfm[3]
    w_o_dead_10hr        : cy.float = cfm[4]
    w_o_dead_100hr       : cy.float = cfm[5]
    w_o_live_herbaceous  : cy.float = cfm[6]
    w_o_live_woody       : cy.float = cfm[7]
    sigma_dead_1hr       : cy.float = cfm[8]
    sigma_dead_10hr      : cy.float = cfm[9]
    sigma_dead_100hr     : cy.float = cfm[10]
    sigma_live_herbaceous: cy.float = cfm[11]
    sigma_live_woody     : cy.float = cfm[12]
    # Expand compressed values
    M_x_dead: cy.float = M_x_dead * 0.01
    h       : cy.float = h * 1000.0
    # Pre-compute some dynamic fuel model values
    dynamic              : cy.bint  = fuel_model_number > 100 and w_o_live_herbaceous > almost_zero
    M_x_dead_herbaceous  : cy.float = M_x_dead              if dynamic else 0.0
    sigma_dead_herbaceous: cy.float = sigma_live_herbaceous if dynamic else 0.0
    # Re-pack everything into a FuelModel struct
    return FuelModel(
        number               = fuel_model_number,
        delta                = delta,
        M_x                  = (M_x_dead,
                                M_x_dead,
                                M_x_dead,
                                M_x_dead_herbaceous,
                                0.0,
                                0.0),
        M_f                  = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
        w_o                  = (w_o_dead_1hr,
                                w_o_dead_10hr,
                                w_o_dead_100hr,
                                0.0,
                                w_o_live_herbaceous,
                                w_o_live_woody),
        sigma                = (sigma_dead_1hr,
                                sigma_dead_10hr,
                                sigma_dead_100hr,
                                sigma_dead_herbaceous,
                                sigma_live_herbaceous,
                                sigma_live_woody),
        h                    = (h, h, h, h, h, h),
        rho_p                = (32.0, 32.0, 32.0, 32.0, 32.0, 32.0),
        S_T                  = (0.0555, 0.0555, 0.0555, 0.0555, 0.0555, 0.0555),
        S_e                  = (0.01, 0.01, 0.01, 0.01, 0.01, 0.01),
        dynamic              = dynamic,
        burnable             = is_burnable_fuel_model_number(fuel_model_number),
        exp_A_sigma          = (compute_exp_A_sigma(-138.0, sigma_dead_1hr),
                                compute_exp_A_sigma(-138.0, sigma_dead_10hr),
                                compute_exp_A_sigma(-138.0, sigma_dead_100hr),
                                compute_exp_A_sigma(-138.0, sigma_dead_herbaceous),
                                compute_exp_A_sigma(-500.0, sigma_live_herbaceous),
                                compute_exp_A_sigma(-500.0, sigma_live_woody)),
        firemod_size_classes = (compute_firemod_size_class(sigma_dead_1hr),
                                compute_firemod_size_class(sigma_dead_10hr),
                                compute_firemod_size_class(sigma_dead_100hr),
                                compute_firemod_size_class(sigma_dead_herbaceous),
                                compute_firemod_size_class(sigma_live_herbaceous),
                                compute_firemod_size_class(sigma_live_woody)),
        f_ij                 = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
        f_i                  = (0.0, 0.0),
        g_ij                 = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
    )


# TODO: OPTIM Replace this dictionary with something more efficient
fuel_model_table = cy.declare(dict, { # dict[int, FuelModel]
    k: expand_compact_fuel_model(k) for k in compact_fuel_model_table.keys()
})


@cy.ccall
@cy.inline
def list_fuel_model_numbers() -> list[int]:
    return list(fuel_model_table.keys())


@cy.ccall
@cy.inline
def list_fuel_models() -> list[FuelModel]:
    return list(fuel_model_table.values())


@cy.ccall
@cy.inline
@cy.exceptval(check=False)
def fuel_model_exists(fuel_model_number: cy.int) -> cy.bint:
    return fuel_model_number in fuel_model_table


@cy.ccall
@cy.inline
@cy.exceptval(check=False)
def get_fuel_model(fuel_model_number: cy.int) -> FuelModel:
    return cy.cast(FuelModel, fuel_model_table[fuel_model_number])

@cy.cfunc
@cy.exceptval(check=False)
def add_dynamic_fuel_loading(fuel_model: FuelModel, M_f: fclaarr) -> FuelModel:
    """
    Updates M_f and w_o. Mutates and returns fuel_model.
    """
    if fuel_model.dynamic:
        # === Dynamic Fuel Model ===
        # Calculate fraction_green and fraction_cured
        w_o                     : fclaarr  = fuel_model.w_o
        live_herbaceous_load    : cy.float = w_o[4]
        live_herbaceous_moisture: cy.float = M_f[4]
        fraction_green          : cy.float = max(0.0, min(1.0, (live_herbaceous_moisture / 0.9) - 0.3333333333333333))
        fraction_cured          : cy.float = 1.0 - fraction_green
        # Set M_f[dead_herbaceous] equal to M_f[dead_1hr]
        M_f_dynamic: fclaarr = (M_f[0], # dead_1hr
                                M_f[1],
                                M_f[2],
                                M_f[0], # dead_herbaceous
                                M_f[4],
                                M_f[5])
        # Set w_o[dead_herbaceous] and w_o[live_herbaceous] based on fraction_cured and fraction_green
        w_o_dynamic: fclaarr = (w_o[0],
                                w_o[1],
                                w_o[2],
                                live_herbaceous_load * fraction_cured, # dead_herbaceous
                                live_herbaceous_load * fraction_green, # live_herbaceous
                                w_o[5])
        # Update the passed in fuel_model to use M_f_dynamic and w_o_dynamic
        fuel_model.M_f = M_f_dynamic
        fuel_model.w_o = w_o_dynamic
        return fuel_model
    else:
        # === Static Fuel Model ===
        # Update the passed in fuel_model to use the passed in M_f
        fuel_model.M_f = M_f
        return fuel_model

# TODO: OPTIM pre-compute this conditional branching since it's fully determined by sigma.
#       This information might be represented efficiently in bit flags.
@cy.cfunc
@cy.exceptval(check=False)
def compute_gij(firemod_size_classes : fclaarr,
                f_ij                 : fclaarr,
                firemod_size_class_ij: cy.float,
                is_dead              : cy.bint) -> cy.float:
    """
    Sums the f_ij of the same category (dead/live) as i, and having the same firemod_size_class.

    NOTE: There may be repetitions in firemod_size_classes, which is why this expression is not
          trivially equal to f_ij[i].
    """
    if is_dead:
        f_ij_0: cy.float = (f_ij[0] if (firemod_size_class_ij == firemod_size_classes[0]) else 0.0)
        f_ij_1: cy.float = (f_ij[1] if (firemod_size_class_ij == firemod_size_classes[1]) else 0.0)
        f_ij_2: cy.float = (f_ij[2] if (firemod_size_class_ij == firemod_size_classes[2]) else 0.0)
        f_ij_3: cy.float = (f_ij[3] if (firemod_size_class_ij == firemod_size_classes[3]) else 0.0)
        return f_ij_0 + f_ij_1 + f_ij_2 + f_ij_3
    else:
        f_ij_4: cy.float = (f_ij[4] if (firemod_size_class_ij == firemod_size_classes[4]) else 0.0)
        f_ij_5: cy.float = (f_ij[5] if (firemod_size_class_ij == firemod_size_classes[5]) else 0.0)
        return f_ij_4 + f_ij_5


@cy.cfunc
@cy.exceptval(check=False)
def add_weighting_factors(fuel_model: FuelModel) -> FuelModel:
    """
    Assigns f_ij, f_i, and g_ij. Mutates and returns fuel_model.
    """
    # Calculate A_ij, A_i, and A_T
    w_o  : fclaarr  = fuel_model.w_o
    sigma: fclaarr  = fuel_model.sigma
    rho_p: fclaarr  = fuel_model.rho_p
    A_ij : fclaarr  = ((sigma[0] * w_o[0]) / rho_p[0],
                       (sigma[1] * w_o[1]) / rho_p[1],
                       (sigma[2] * w_o[2]) / rho_p[2],
                       (sigma[3] * w_o[3]) / rho_p[3],
                       (sigma[4] * w_o[4]) / rho_p[4],
                       (sigma[5] * w_o[5]) / rho_p[5]) # TODO: OPTIM pre-compute sigma/rho_p
    A_i_0: cy.float = A_ij[0] + A_ij[1] + A_ij[2] + A_ij[3]
    A_i_1: cy.float = A_ij[4] + A_ij[5]
    A_i  : fcatarr  = (A_i_0, A_i_1)
    A_T  : cy.float = A_i_0 + A_i_1
    # Calculate f_ij
    f_ij_0   : cy.float = 0.0
    f_ij_1   : cy.float = 0.0
    f_ij_2   : cy.float = 0.0
    f_ij_3   : cy.float = 0.0
    f_ij_4   : cy.float = 0.0
    f_ij_5   : cy.float = 0.0
    A_i_0_inv: cy.float
    A_i_1_inv: cy.float
    if A_i_0 > almost_zero:
        A_i_0_inv = 1.0 / A_i_0
        f_ij_0    = A_ij[0] * A_i_0_inv
        f_ij_1    = A_ij[1] * A_i_0_inv
        f_ij_2    = A_ij[2] * A_i_0_inv
        f_ij_3    = A_ij[3] * A_i_0_inv
    if A_i_1 > almost_zero:
        A_i_1_inv = 1.0 / A_i_1
        f_ij_4    = A_ij[4] * A_i_1_inv
        f_ij_5    = A_ij[5] * A_i_1_inv
    f_ij: fclaarr = (f_ij_0, f_ij_1, f_ij_2, f_ij_3, f_ij_4, f_ij_5)
    # Calculate f_i
    f_i_0  : cy.float = 0.0
    f_i_1  : cy.float = 0.0
    A_T_inv: cy.float
    if A_T > almost_zero:
        A_T_inv = 1.0 / A_T
        f_i_0   = A_i_0 * A_T_inv
        f_i_1   = A_i_1 * A_T_inv
    f_i: fcatarr = (f_i_0, f_i_1)
    # Calculate g_ij
    firemod_size_classes: fclaarr = fuel_model.firemod_size_classes
    g_ij                : fclaarr = (compute_gij(firemod_size_classes, f_ij, firemod_size_classes[0], True),
                                     compute_gij(firemod_size_classes, f_ij, firemod_size_classes[1], True),
                                     compute_gij(firemod_size_classes, f_ij, firemod_size_classes[2], True),
                                     compute_gij(firemod_size_classes, f_ij, firemod_size_classes[3], True),
                                     compute_gij(firemod_size_classes, f_ij, firemod_size_classes[4], False),
                                     compute_gij(firemod_size_classes, f_ij, firemod_size_classes[5], False))
    # Update the passed in fuel_model to use f_ij, f_i, and g_ij
    fuel_model.f_ij = f_ij
    fuel_model.f_i  = f_i
    fuel_model.g_ij = g_ij
    return fuel_model

@cy.cfunc
@cy.exceptval(check=False)
def add_live_moisture_of_extinction(fuel_model: FuelModel) -> FuelModel:
    """
    Equation 88 from Rothermel 1972 adjusted by Albini 1976 Appendix III.

    Updates M_x. Mutates and returns fuel_model.
    """
    # Calculate dead_moisture_factor, dead_loading_factor, and live_loading_factor
    w_o                 : fclaarr  = fuel_model.w_o
    exp_A_sigma         : fclaarr  = fuel_model.exp_A_sigma
    M_f                 : fclaarr  = fuel_model.M_f
    M_x                 : fclaarr  = fuel_model.M_x
    M_x_dead            : cy.float = M_x[0]
    loading_factors     : fclaarr  = (w_o[0] * exp_A_sigma[0],
                                      w_o[1] * exp_A_sigma[1],
                                      w_o[2] * exp_A_sigma[2],
                                      w_o[3] * exp_A_sigma[3],
                                      w_o[4] * exp_A_sigma[4],
                                      w_o[5] * exp_A_sigma[5])
    dead_moisture_factor: cy.float = (M_f[0] * loading_factors[0] +
                                      M_f[1] * loading_factors[1] +
                                      M_f[2] * loading_factors[2] +
                                      M_f[3] * loading_factors[3])
    dead_loading_factor : cy.float = loading_factors[0] + loading_factors[1] + loading_factors[2] + loading_factors[3]
    live_loading_factor : cy.float = loading_factors[4] + loading_factors[5]
    # Calculate M_x_live
    dead_fuel_moisture: cy.float
    dead_to_live_ratio: cy.float
    M_x_live          : cy.float
    if (dead_loading_factor > almost_zero and live_loading_factor > almost_zero):
        dead_fuel_moisture = dead_moisture_factor / dead_loading_factor
        dead_to_live_ratio = dead_loading_factor / live_loading_factor
        M_x_live           = max(M_x_dead,
                                 (2.9 * dead_to_live_ratio * (1.0 - (dead_fuel_moisture / M_x_dead))) - 0.226)
    else:
        M_x_live = M_x_dead
    # Calculate M_x_new
    M_x_new: fclaarr = (M_x[0],
                        M_x[1],
                        M_x[2],
                        M_x[3],
                        M_x_live,
                        M_x_live)
    # Update the passed in fuel_model to use M_x_new
    fuel_model.M_x = M_x_new
    return fuel_model

@cy.ccall
@cy.exceptval(check=False)
def moisturize(fuel_model: FuelModel, fuel_moisture: fclaarr) -> FuelModel:
    """
    Updates w_o, M_f, and M_x and assigns f_ij, f_i, and g_ij.
    Returns a new FuelModel struct.
    """
    # Create a copy of fuel_model
    fuel_model_copy: FuelModel = FuelModel(
        number               = fuel_model.number,
        delta                = fuel_model.delta,
        M_x                  = fuel_model.M_x,
        M_f                  = fuel_model.M_f,
        w_o                  = fuel_model.w_o,
        sigma                = fuel_model.sigma,
        h                    = fuel_model.h,
        rho_p                = fuel_model.rho_p,
        S_T                  = fuel_model.S_T,
        S_e                  = fuel_model.S_e,
        dynamic              = fuel_model.dynamic,
        burnable             = fuel_model.burnable,
        exp_A_sigma          = fuel_model.exp_A_sigma,
        firemod_size_classes = fuel_model.firemod_size_classes,
        f_ij                 = fuel_model.f_ij,
        f_i                  = fuel_model.f_i,
        g_ij                 = fuel_model.g_ij,
    )
    # Mutate and return the copy
    dynamic_fuel_model    : FuelModel = add_dynamic_fuel_loading(fuel_model_copy, fuel_moisture)
    weighted_fuel_model   : FuelModel = add_weighting_factors(dynamic_fuel_model)
    moisturized_fuel_model: FuelModel = add_live_moisture_of_extinction(weighted_fuel_model)
    return moisturized_fuel_model

from pprint import pprint
from pyretechnics.fuel_models import get_fuel_model

static_fuel_model_number  = 1   # R01
dynamic_fuel_model_number = 101 # GR1

static_fuel_model  = get_fuel_model(static_fuel_model_number)
dynamic_fuel_model = get_fuel_model(dynamic_fuel_model_number)

print("Static Fuel Model: R01")
pprint(static_fuel_model)

print("\nDynamic Fuel Model: GR1")
pprint(dynamic_fuel_model)

from pprint import pprint
from pyretechnics.fuel_models import get_fuel_model, moisturize

static_fuel_model_number  = 1   # R01
dynamic_fuel_model_number = 101 # GR1

static_fuel_model  = get_fuel_model(static_fuel_model_number)
dynamic_fuel_model = get_fuel_model(dynamic_fuel_model_number)

fuel_moisture = (
    0.05, # dead_1hr
    0.10, # dead_10hr
    0.15, # dead_100hr
    0.00, # dead_herbaceous (will be set by moisturize for dynamic fuel models)
    0.30, # live_herbaceous
    0.50, # live_woody
)

print("Moisturized Static Fuel Model: R01")
pprint(moisturize(static_fuel_model, fuel_moisture))

print("\nMoisturized Dynamic Fuel Model: GR1")
pprint(moisturize(dynamic_fuel_model, fuel_moisture))

import cython
import cython as cy
if cython.compiled:
    from cython.cimports.libc.math import sqrt, pow, log, exp
    from cython.cimports.pyretechnics.cy_types import \
        vec_xy, vec_xyz, fcatarr, fclaarr, FuelModel, ProjectedVectors, FireBehaviorMin, \
        FireBehaviorMax, SpreadBehavior
    import cython.cimports.pyretechnics.conversion as conv
    import cython.cimports.pyretechnics.vector_utils as vu
else:
    from math import sqrt, pow, log, exp
    from pyretechnics.py_types import \
        vec_xy, vec_xyz, fcatarr, fclaarr, FuelModel, ProjectedVectors, FireBehaviorMin, \
        FireBehaviorMax, SpreadBehavior
    import pyretechnics.conversion as conv
    import pyretechnics.vector_utils as vu

@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def calc_residence_time(sigma_prime: cy.float) -> cy.float:
    """
    Returns the residence time (total burning time) of fuel (min) given:
    - sigma_prime :: ft^2/ft^3 (surface area to volume ratio)
    """
    return 384.0 / sigma_prime if (sigma_prime > 0.0) else 0.0


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def calc_flame_depth(spread_rate: cy.float, residence_time: cy.float) -> cy.float:
    """
    Returns the depth, or front-to-back distance, of the actively flaming zone
    of a free-spreading fire (ft) given:
    - spread_rate    :: ft/min (orthogonal to the fireline)
    - residence_time :: min
    """
    return spread_rate * residence_time


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def calc_fireline_intensity(reaction_intensity: cy.float, flame_depth: cy.float) -> cy.float:
    """
    Returns the rate of heat release per unit of fire edge (Btu/ft/s) given:
    - reaction_intensity :: Btu/ft^2/min
    - flame_depth        :: ft
    """
    return (reaction_intensity * flame_depth) / 60.0


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def calc_flame_length(fireline_intensity: cy.float) -> cy.float:
    """
    Returns the average flame length (m) given:
    - fireline_intensity :: kW/m
    """
    return 0.07747042253266703 * pow(fireline_intensity, 0.46)


@cy.ccall
@cy.inline
@cy.exceptval(check=False)
def calc_areal_heat_output(spread_rate: cy.float, fireline_intensity: cy.float) -> cy.float:
    """
    Returns the heat per unit area (kJ/m^2) given:
    - spread_rate        :: m/min
    - fireline_intensity :: kW/m
    """
    return 60.0 * fireline_intensity / spread_rate if spread_rate > 0.0 else 0.0

@cy.cfunc
@cy.exceptval(check=False)
def make_surface_fire_min(base_spread_rate        : cy.float,
                          base_fireline_intensity : cy.float,
                          max_effective_wind_speed: cy.float,
                          B                       : cy.float,
                          C                       : cy.float,
                          F                       : cy.float,
                          beta                    : cy.float) -> FireBehaviorMin:
    _phiS_G: cy.float = 0.0
    if (beta > 0.0):
        _phiS_G = 5.275 * pow(beta, -0.3)

    _phiW_scalr: cy.float = 0.0
    _phiW_expnt: cy.float = 0.0
    if (F > 0.0):
        _phiW_scalr = (C / F) * pow(conv.m_to_ft(1.0), B)
        _phiW_expnt = B

    _ws_scalr: cy.float = 0.0
    _ws_expnt: cy.float = 0.0
    if (B > 0.0):
        B_inverse: cy.float = 1.0 / B
        _ws_scalr = conv.ft_to_m(pow((F / C), B_inverse))
        _ws_expnt = B_inverse

    return FireBehaviorMin(
        base_spread_rate,
        base_fireline_intensity,
        max_effective_wind_speed,
        _phiS_G,
        _phiW_scalr,
        _phiW_expnt,
        _ws_scalr,
        _ws_expnt)


@cy.ccall
@cy.exceptval(check=False)
def calc_surface_fire_behavior_no_wind_no_slope(moisturized_fuel_model: FuelModel,
                                                spread_rate_adjustment: cy.float = 1.0) -> FireBehaviorMin:
    """
    Given these inputs:
    - moisturized_fuel_model :: FuelModel struct of fuel model and fuel moisture properties
      - delta :: ft                                        (fuel depth)
      - w_o   :: lb/ft^2                                   (ovendry fuel loading)
      - rho_p :: lb/ft^3                                   (ovendry particle density)
      - sigma :: ft^2/ft^3                                 (fuel particle surface-area-to-volume ratio)
      - h     :: Btu/lb                                    (fuel particle low heat content)
      - S_T   :: lb minerals/lb ovendry weight             (fuel particle total mineral content)
      - S_e   :: lb silica-free minerals/lb ovendry weight (fuel particle effective mineral content)
      - M_x   :: lb moisture/lb ovendry weight             (fuel particle moisture of extinction)
      - M_f   :: lb moisture/lb ovendry weight             (fuel particle moisture content)
      - f_ij  :: %                                         (percent load per size class)
      - f_i   :: %                                         (percent load per category)
      - g_ij  :: %                                         (percent load per size class - Albini_1976_FIREMOD, page 20)
    - spread_rate_adjustment :: unitless float (1.0 for no adjustment)

    return a FireBehaviorMin struct containing these keys:
    - base_spread_rate         :: m/min
    - base_fireline_intensity  :: kW/m
    - max_effective_wind_speed :: m/min
    - _phiS_G                  :: intermediate value for computing phi_S (unitless)
    - _phiW_scalr              :: intermediate value for computing phi_W (unitless)
    - _phiW_expnt              :: intermediate value for computing phi_W (unitless)
    - _ws_scalr                :: intermediate value for computing effective_wind_speed (m/min)
    - _ws_expnt                :: intermediate value for computing effective_wind_speed (unitless)
    """
    # Unpack fuel model values
    delta: cy.float = moisturized_fuel_model.delta
    w_o  : fclaarr  = moisturized_fuel_model.w_o
    rho_p: fclaarr  = moisturized_fuel_model.rho_p
    sigma: fclaarr  = moisturized_fuel_model.sigma
    M_f  : fclaarr  = moisturized_fuel_model.M_f
    f_ij : fclaarr  = moisturized_fuel_model.f_ij
    f_i  : fcatarr  = moisturized_fuel_model.f_i
    # Calculate base spread rate (no wind, no slope)
    sigma_prime: cy.float = calc_surface_area_to_volume_ratio(f_i, f_ij, sigma)
    beta       : cy.float = calc_packing_ratio(w_o, rho_p, delta)
    beta_op    : cy.float = calc_optimum_packing_ratio(sigma_prime)
    I_R        : cy.float = calc_reaction_intensity(moisturized_fuel_model, sigma_prime, beta, beta_op) # Btu/ft^2/min
    xi         : cy.float = calc_propagating_flux_ratio(sigma_prime, beta)
    heat_source: cy.float = calc_heat_source(I_R, xi)                             # Btu/ft^2/min
    rho_b      : cy.float = calc_ovendry_bulk_density(w_o, delta)                 # lb/ft^3
    epsilon_ij : fclaarr  = calc_effective_heating_number_distribution(sigma)
    Q_ig_ij    : fclaarr  = calc_heat_of_preignition_distribution(M_f)            # Btu/lb
    heat_sink  : cy.float = calc_heat_sink(f_i, f_ij, rho_b, epsilon_ij, Q_ig_ij) # Btu/ft^3
    R0         : cy.float = calc_spread_rate(heat_source, heat_sink)              # ft/min
    # Calculate base fireline intensity (no wind, no slope)
    t_res: cy.float = calc_residence_time(sigma_prime)  # min
    D_A  : cy.float = calc_flame_depth(R0, t_res)       # ft
    I_s  : cy.float = calc_fireline_intensity(I_R, D_A) # Btu/ft/s
    # Pre-compute values related to wind and slope
    U_eff_max: cy.float = calc_max_effective_wind_speed(I_R) # ft/min
    B        : cy.float = 0.02526 * pow(sigma_prime, 0.54)
    C        : cy.float = 7.47 * exp(-0.133 * pow(sigma_prime, 0.55))
    E        : cy.float = 0.715 * exp(-3.59 * (sigma_prime * 1e-4))
    F        : cy.float = pow((beta / beta_op), E)
    # Return no-wind-no-slope surface fire behavior values
    base_spread_rate        : cy.float = conv.ft_to_m(R0 * spread_rate_adjustment)
    base_fireline_intensity : cy.float = conv.Btu_ft_s_to_kW_m(I_s * spread_rate_adjustment)
    max_effective_wind_speed: cy.float = conv.ft_to_m(U_eff_max)
    return make_surface_fire_min(
        base_spread_rate,
        base_fireline_intensity,
        max_effective_wind_speed,
        B,
        C,
        F,
        beta
    )

@cy.cfunc
@cy.exceptval(check=False)
def calc_wind_adjustment_factor(fuel_bed_depth: cy.float, canopy_height: cy.float, canopy_cover: cy.float) -> cy.float:
    """
    Return the wind adjustment factor (unitless) given these inputs:
    - fuel_bed_depth :: ft
    - canopy_height  :: ft
    - canopy_cover   :: 0-1
    """
    if (canopy_cover > 0.0) and (canopy_height > 0.0):
        # sheltered: equation 2 based on CC and CH, CR=1 (Andrews 2012)
        A: cy.float = sqrt((canopy_cover / 3.0) * canopy_height)
        B: cy.float = log((20.0 + 0.36 * canopy_height) / (0.13 * canopy_height))
        return 0.555 / (A * B)
    elif (fuel_bed_depth > 0.0):
        # unsheltered: equation 6 H_F = H (Andrews 2012)
        A: cy.float = log((20.0 + 0.36 * fuel_bed_depth) / (0.13 * fuel_bed_depth))
        return 1.83 / A # 1.83 truncated from 1.8328795184533409
    else:
        # non-burnable fuel model
        return 0.0


@cy.ccall
@cy.exceptval(check=False)
def calc_midflame_wind_speed(wind_speed_20ft: cy.float,
                             fuel_bed_depth : cy.float,
                             canopy_height  : cy.float,
                             canopy_cover   : cy.float) -> cy.float:
    """
    Return the midflame wind speed (S) given these inputs:
    - wind_speed_20ft :: S
    - fuel_bed_depth  :: ft
    - canopy_height   :: ft
    - canopy_cover    :: 0-1
    """
    wind_adj_factor: cy.float = calc_wind_adjustment_factor(fuel_bed_depth, canopy_height, canopy_cover)
    return wind_speed_20ft * wind_adj_factor

@cy.cfunc
@cy.exceptval(check=False)
def project_wind_and_slope_vectors_3d(wind_speed        : cy.float,
                                      downwind_direction: cy.float,
                                      slope             : cy.float,
                                      upslope_direction : cy.float) -> ProjectedVectors:
    """
    Given these inputs:
    - wind_speed         :: S
    - downwind_direction :: degrees clockwise from North
    - slope              :: rise/run
    - upslope_direction  :: degrees clockwise from North

    return a ProjectedVectors struct containing these keys:
    - wind_vector_3d  :: (x: S, y: S, z: S)
    - slope_vector_3d :: (x, y, z)
    """
    # Convert wind and slope vectors from azimuthal to cartesian coordinates
    wind_vector_2d : vec_xy = conv.azimuthal_to_cartesian(wind_speed, downwind_direction)
    slope_vector_2d: vec_xy = conv.azimuthal_to_cartesian(slope, upslope_direction)
    # Project wind and slope vectors onto the slope-tangential plane
    wind_vector_3d : vec_xyz = vu.to_slope_plane(wind_vector_2d, slope_vector_2d)
    slope_vector_3d: vec_xyz = vu.to_slope_plane(slope_vector_2d, slope_vector_2d)
    return ProjectedVectors(wind_vector_3d, slope_vector_3d)


@cy.cfunc
@cy.exceptval(check=False)
def get_phi_E(wind_vector_3d: vec_xyz, slope_vector_3d: vec_xyz, phi_W: cy.float, phi_S: cy.float) -> vec_xyz:
    # Convert wind and slope vectors to unit vectors on the slope-tangential plane
    w_S: vec_xyz = vu.as_unit_vector_3d(wind_vector_3d)  if phi_W > 0.0 else wind_vector_3d
    u_S: vec_xyz = vu.as_unit_vector_3d(slope_vector_3d) if phi_S > 0.0 else slope_vector_3d
    # Create the 3D slope-tangential phi_W, phi_S, and phi_E vectors
    phi_W_3d: vec_xyz = vu.scale_3d(phi_W, w_S)
    phi_S_3d: vec_xyz = vu.scale_3d(phi_S, u_S)
    phi_E_3d: vec_xyz = vu.add_3d(phi_W_3d, phi_S_3d)
    return phi_E_3d

# TODO: Change model from str to enumerated type
@cy.cfunc
def surface_length_to_width_ratio(effective_wind_speed: cy.float, model: str = "behave") -> cy.float:
    """
    Calculate the length_to_width_ratio of the surface fire front given:
    - effective_wind_speed :: m/min (aligned with the slope-tangential plane)
    - model                :: "rothermel" or "behave" (Optional)
    """
    effective_wind_speed_mph = conv.m_min_to_mph(effective_wind_speed)
    # Select formula by model
    if model == "rothermel":
        return 1.0 + 0.25 * effective_wind_speed_mph

    elif model == "behave":
        return min(8.0,
                   0.936 * exp(0.1147 * effective_wind_speed_mph)
                   +
                   0.461 * exp(-0.0692 * effective_wind_speed_mph)
                   -
                   0.397)

    else:
        raise ValueError("Invalid input: model must be 'rothermel' or 'behave'.")


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def surface_fire_eccentricity(length_to_width_ratio: cy.float) -> cy.float:
    """
    Calculate the eccentricity (E) of the surface fire front using eq. 8 from
    Albini and Chase 1980 given:
    - L/W :: (1: circular spread, > 1: elliptical spread)
    """
    return sqrt((length_to_width_ratio * length_to_width_ratio) - 1.0) / length_to_width_ratio

@cy.cfunc
@cy.exceptval(check=False)
def maybe_limit_wind_speed(use_wind_limit  : cy.bint,
                           max_wind_speed  : cy.float,
                           surface_fire_min: FireBehaviorMin,
                           phi_E_magnitude : cy.float) -> tuple[cy.float, cy.float]:
    """
    Given these inputs:
    - use_wind_limit   :: boolean
    - max_wind_speed   :: m/min
    - surface_fire_min :: FireBehaviorMin struct of no-wind-no-slope surface fire behavior values
      - base_spread_rate         :: m/min
      - base_fireline_intensity  :: kW/m
      - max_effective_wind_speed :: m/min
      - _phiS_G                  :: intermediate value for computing phi_S (unitless)
      - _phiW_scalr              :: intermediate value for computing phi_W (unitless)
      - _phiW_expnt              :: intermediate value for computing phi_W (unitless)
      - _ws_scalr                :: intermediate value for computing effective_wind_speed (m/min)
      - _ws_expnt                :: intermediate value for computing effective_wind_speed (unitless)
    - phi_E_magnitude  :: unitless

    return a tuple with these fields:
    - limited_wind_speed :: m/min
    - limited_phi_E      :: unitless
    """
    effective_wind_speed: cy.float = get_wind_speed(surface_fire_min, phi_E_magnitude)
    if (use_wind_limit and effective_wind_speed > max_wind_speed):
        return (
            max_wind_speed,
            get_phi_W(surface_fire_min, max_wind_speed),
        )
    else:
        return (
            effective_wind_speed,
            phi_E_magnitude,
        )


@cy.ccall
@cy.exceptval(check=False)
def calc_surface_fire_behavior_max(surface_fire_min      : FireBehaviorMin,
                                   midflame_wind_speed   : cy.float,
                                   upwind_direction      : cy.float,
                                   slope                 : cy.float,
                                   aspect                : cy.float,
                                   use_wind_limit        : cy.bint = True,
                                   surface_lw_ratio_model: str = "behave") -> FireBehaviorMax:
    """
    Given these inputs:
    - surface_fire_min       :: FireBehaviorMin struct of no-wind-no-slope surface fire behavior values
      - base_spread_rate         :: m/min
      - base_fireline_intensity  :: kW/m
      - max_effective_wind_speed :: m/min
      - _phiS_G                  :: intermediate value for computing phi_S (unitless)
      - _phiW_scalr              :: intermediate value for computing phi_W (unitless)
      - _phiW_expnt              :: intermediate value for computing phi_W (unitless)
      - _ws_scalr                :: intermediate value for computing effective_wind_speed (m/min)
      - _ws_expnt                :: intermediate value for computing effective_wind_speed (unitless)
    - midflame_wind_speed    :: m/min
    - upwind_direction       :: degrees clockwise from North
    - slope                  :: rise/run
    - aspect                 :: degrees clockwise from North
    - use_wind_limit         :: boolean (Optional)
    - surface_lw_ratio_model :: "rothermel" or "behave" (Optional)

    return a FireBehaviorMax struct containing these keys:
    - max_fire_type          :: 1 (surface)
    - max_spread_rate        :: m/min
    - max_spread_direction   :: (x, y, z) unit vector
    - max_fireline_intensity :: kW/m
    - max_flame_length       :: m
    - length_to_width_ratio  :: unitless (1: circular spread, > 1: elliptical spread)
    - eccentricity           :: unitless (0: circular spread, > 0: elliptical spread)
    - critical_spread_rate   :: m/min
    """
    # Unpack no-wind-no-slope surface fire behavior values
    spread_rate       : cy.float = surface_fire_min.base_spread_rate
    fireline_intensity: cy.float = surface_fire_min.base_fireline_intensity
    max_wind_speed    : cy.float = surface_fire_min.max_effective_wind_speed
    # Reverse the provided wind and slope directions
    downwind_direction: cy.float = conv.opposite_direction(upwind_direction)
    upslope_direction : cy.float = conv.opposite_direction(aspect)
    # Project wind and slope vectors onto the slope-tangential plane
    vectors        : ProjectedVectors = project_wind_and_slope_vectors_3d(midflame_wind_speed,
                                                                          downwind_direction,
                                                                          slope,
                                                                          upslope_direction)
    wind_vector_3d : vec_xyz = vectors.wind_vector_3d  # m/min
    slope_vector_3d: vec_xyz = vectors.slope_vector_3d # rise/run
    # Calculate phi_W and phi_S
    # NOTE: |wind_vector_3d| = slope-aligned midflame wind speed
    phi_W: cy.float = get_phi_W(surface_fire_min, vu.vector_magnitude_3d(wind_vector_3d))
    phi_S: cy.float = get_phi_S(surface_fire_min, slope)
    # Calculate phi_E and the max_spread_direction
    phi_E_3d: vec_xyz  = get_phi_E(wind_vector_3d, slope_vector_3d, phi_W, phi_S)
    phi_E   : cy.float = vu.vector_magnitude_3d(phi_E_3d)
    max_spread_direction: vec_xyz
    if phi_E > 0.0:
        max_spread_direction = vu.as_unit_vector_3d(phi_E_3d)
    elif phi_S > 0.0:
        max_spread_direction = vu.as_unit_vector_3d(slope_vector_3d)
    else:
        max_spread_direction = (0.0, 1.0, 0.0) # default: North
    # Limit effective wind speed to max wind speed if use_wind_limit == True
    (limited_wind_speed, limited_phi_E) = maybe_limit_wind_speed(use_wind_limit,
                                                                 max_wind_speed,
                                                                 surface_fire_min,
                                                                 phi_E)
    # Calculate and return max surface fire behavior values
    max_spread_rate       : cy.float = spread_rate * (1.0 + limited_phi_E)
    max_fireline_intensity: cy.float = fireline_intensity * (1.0 + limited_phi_E)
    length_to_width_ratio : cy.float = surface_length_to_width_ratio(limited_wind_speed, surface_lw_ratio_model)
    return FireBehaviorMax(
        max_fire_type          = 1,
        max_spread_rate        = max_spread_rate,
        max_spread_direction   = max_spread_direction,
        max_fireline_intensity = max_fireline_intensity,
        max_flame_length       = calc_flame_length(max_fireline_intensity),
        length_to_width_ratio  = length_to_width_ratio,
        eccentricity           = surface_fire_eccentricity(length_to_width_ratio),
        critical_spread_rate   = 0.0,
    )

@cy.ccall
@cy.exceptval(check=False)
def calc_surface_fire_behavior_in_direction(surface_fire_max: FireBehaviorMax,
                                            spread_direction: vec_xyz) -> SpreadBehavior:
    """
    Given these inputs:
    - surface_fire_max     :: FireBehaviorMax struct of max surface fire behavior values
      - max_fire_type          :: 1 (surface)
      - max_spread_rate        :: m/min
      - max_spread_direction   :: (x, y, z) unit vector
      - max_fireline_intensity :: kW/m
      - max_flame_length       :: m
      - length_to_width_ratio  :: unitless (1: circular spread, > 1: elliptical spread)
      - eccentricity           :: unitless (0: circular spread, > 0: elliptical spread)
      - critical_spread_rate   :: m/min
    - spread_direction     :: 3D unit vector on the slope-tangential plane

    return a SpreadBehavior struct containing these keys:
    - dphi_dt            :: phi/min
    - fire_type          :: 1 (surface)
    - spread_rate        :: m/min
    - spread_direction   :: (x, y, z) unit vector
    - fireline_intensity :: kW/m
    - flame_length       :: m
    """
    # Unpack max surface fire behavior values
    max_spread_rate       : cy.float = surface_fire_max.max_spread_rate
    max_spread_direction  : vec_xyz  = surface_fire_max.max_spread_direction
    max_fireline_intensity: cy.float = surface_fire_max.max_fireline_intensity
    eccentricity          : cy.float = surface_fire_max.eccentricity
    # Calculate cos(w), where w is the offset angle between these unit vectors on the slope-tangential plane
    cos_w: cy.float = vu.dot_3d(max_spread_direction, spread_direction)
    # Calculate adjustment due to the offset angle from the max spread direction
    adjustment: cy.float = (1.0 - eccentricity) / (1.0 - eccentricity * cos_w)
    # Update surface fire behavior values by the adjustment value
    fireline_intensity: cy.float = max_fireline_intensity * adjustment
    return SpreadBehavior(
        dphi_dt            = 0.0,
        fire_type          = 1, # surface
        spread_rate        = max_spread_rate * adjustment,
        spread_direction   = spread_direction,
        fireline_intensity = fireline_intensity,
        flame_length       = calc_flame_length(fireline_intensity),
    )

from pprint import pprint
import pyretechnics.fuel_models as fm
import pyretechnics.surface_fire as sf

# Set input parameters
fuel_model    = fm.get_fuel_model(101) # GR1
fuel_moisture = (
    0.05, # dead_1hr
    0.10, # dead_10hr
    0.15, # dead_100hr
    0.00, # dead_herbaceous (will be set by moisturize for dynamic fuel models)
    0.90, # live_herbaceous
    0.60, # live_woody
)

# Apply fuel moisture to fuel model
moisturized_fuel_model = fm.moisturize(fuel_model, fuel_moisture)

# Calculate no-wind-no-slope surface fire behavior
surface_fire_min = sf.calc_surface_fire_behavior_no_wind_no_slope(moisturized_fuel_model)

# View results
print("No-Wind-No-Slope Surface Fire Behavior for Fuel Model GR1 with Fuel Moisture = "
      + str(moisturized_fuel_model["M_f"]))
pprint({k: surface_fire_min[k] for k in ["base_spread_rate", "base_fireline_intensity", "max_effective_wind_speed"]})

import pyretechnics.conversion as conv
import pyretechnics.fuel_models as fm
import pyretechnics.surface_fire as sf

# Set input parameters
fuel_model       = fm.get_fuel_model(101) # GR1
canopy_height    = 20                     # meters
canopy_cover     = 0.6                    # 0-1
wind_speed_10m   = 20                     # km/hr
upwind_direction = 45.0                   # degrees clockwise from North

# Convert from 10m wind speed to 20ft wind speed
wind_speed_20ft = conv.wind_speed_10m_to_wind_speed_20ft(wind_speed_10m) # km/hr

# Convert from 20ft wind speed to midflame wind speed in m/min
midflame_wind_speed = sf.calc_midflame_wind_speed(conv.km_hr_to_m_min(wind_speed_20ft), # m/min
                                                  fuel_model["delta"],                  # ft
                                                  conv.m_to_ft(canopy_height),          # ft
                                                  canopy_cover)                         # 0-1

# View results
print("Midflame Wind Speed and Direction for 10m Wind Speed = 20 km/hr, Upwind Direction = 45.0\n"
      + " with Fuel Model GR1, Canopy Height = 20m, and Canopy Cover = 60%\n"
      + "- Speed: " + str(midflame_wind_speed) + " (m/min)\n"
      + "- Downwind Direction: " + str(conv.opposite_direction(upwind_direction)) + " (degrees clockwise from North)")

from pprint import pprint
import pyretechnics.fuel_models as fm
import pyretechnics.surface_fire as sf

# Set input parameters
fuel_model    = fm.get_fuel_model(101) # GR1
fuel_moisture = (
    0.05, # dead_1hr
    0.10, # dead_10hr
    0.15, # dead_100hr
    0.00, # dead_herbaceous (will be set by moisturize for dynamic fuel models)
    0.90, # live_herbaceous
    0.60, # live_woody
)
midflame_wind_speed = 500.0 # m/min
upwind_direction    = 215.0 # degrees clockwise from North
slope               = 0.2   # rise/run
aspect              = 270.0 # degrees clockwise from North

# Apply fuel moisture to fuel model
moisturized_fuel_model = fm.moisturize(fuel_model, fuel_moisture)

# Calculate no-wind-no-slope surface fire behavior
surface_fire_min = sf.calc_surface_fire_behavior_no_wind_no_slope(moisturized_fuel_model)

# Calculate surface fire behavior in the direction of maximum spread with effective wind speed limit
surface_fire_max = sf.calc_surface_fire_behavior_max(surface_fire_min,
                                                     midflame_wind_speed,
                                                     upwind_direction,
                                                     slope,
                                                     aspect,
                                                     surface_lw_ratio_model="rothermel")

# View results
print("Surface Fire Behavior Max (Limited) for Fuel Model GR1 with Fuel Moisture = "
      + str(moisturized_fuel_model["M_f"]))
pprint(surface_fire_max)

# Calculate surface fire behavior in the direction of maximum spread without effective wind speed limit
surface_fire_max_unlimited = sf.calc_surface_fire_behavior_max(surface_fire_min,
                                                               midflame_wind_speed,
                                                               upwind_direction,
                                                               slope,
                                                               aspect,
                                                               use_wind_limit=False,
                                                               surface_lw_ratio_model="rothermel")

# View results
print("\nSurface Fire Behavior Max (Unlimited) for Fuel Model GR1 with Fuel Moisture = "
      + str(moisturized_fuel_model["M_f"]))
pprint(surface_fire_max_unlimited)

import numpy as np
from pprint import pprint
import pyretechnics.fuel_models as fm
import pyretechnics.surface_fire as sf
import pyretechnics.vector_utils as vu

# Set input parameters
fuel_model    = fm.get_fuel_model(101) # GR1
fuel_moisture = (
    0.05, # dead_1hr
    0.10, # dead_10hr
    0.15, # dead_100hr
    0.00, # dead_herbaceous (will be set by moisturize for dynamic fuel models)
    0.90, # live_herbaceous
    0.60, # live_woody
)
midflame_wind_speed = 500.0 # m/min
upwind_direction    = 215.0 # degrees clockwise from North
slope               = 0.2   # rise/run
aspect              = 270.0 # degrees clockwise from North

# Apply fuel moisture to fuel model
moisturized_fuel_model = fm.moisturize(fuel_model, fuel_moisture)

# Calculate no-wind-no-slope surface fire behavior
surface_fire_min = sf.calc_surface_fire_behavior_no_wind_no_slope(moisturized_fuel_model)

# Calculate surface fire behavior in the direction of maximum spread with effective wind speed limit
surface_fire_max = sf.calc_surface_fire_behavior_max(surface_fire_min,
                                                     midflame_wind_speed,
                                                     upwind_direction,
                                                     slope,
                                                     aspect,
                                                     surface_lw_ratio_model="rothermel")

# Calculate surface fire behavior for the heading, flanking, and backing fires
heading_fire_direction  = surface_fire_max["max_spread_direction"]
# FIXME: Change upstream functions to return max_spread_direction as type vec_xyz
#        and change downstream functions to accept spread_direction as type vec_xyz
(dx, dy, dz)            = heading_fire_direction
flanking_fire_direction = np.asarray(vu.rotate_on_sloped_plane((dx, dy, dz), 90.0, slope, aspect))
backing_fire_direction  = np.asarray(vu.rotate_on_sloped_plane((dx, dy, dz), 180.0, slope, aspect))

heading_fire_behavior  = sf.calc_surface_fire_behavior_in_direction(surface_fire_max, heading_fire_direction)
flanking_fire_behavior = sf.calc_surface_fire_behavior_in_direction(surface_fire_max, flanking_fire_direction)
backing_fire_behavior  = sf.calc_surface_fire_behavior_in_direction(surface_fire_max, backing_fire_direction)

# View heading fire results
print("Heading Fire Behavior for Fuel Model GR1 with Fuel Moisture = "
      + str(moisturized_fuel_model["M_f"]))
pprint(heading_fire_behavior)

# View flanking fire results
print("\nFlanking Fire Behavior for Fuel Model GR1 with Fuel Moisture = "
      + str(moisturized_fuel_model["M_f"]))
pprint(flanking_fire_behavior)

# View backing fire results
print("\nBacking Fire Behavior for Fuel Model GR1 with Fuel Moisture = "
      + str(moisturized_fuel_model["M_f"]))
pprint(backing_fire_behavior)

import cython
import cython as cy
if cython.compiled:
    from cython.cimports.libc.math import sqrt, exp, pow
    from cython.cimports.pyretechnics.cy_types import \
        vec_xyz, ProjectedVectors, FireBehaviorMax, SpreadBehavior, CrownSpreadInfo
    import cython.cimports.pyretechnics.conversion as conv
    import cython.cimports.pyretechnics.vector_utils as vu
    import cython.cimports.pyretechnics.surface_fire as sf
else:
    from math import sqrt, exp, pow
    from pyretechnics.py_types import \
        vec_xyz, ProjectedVectors, FireBehaviorMax, SpreadBehavior, CrownSpreadInfo
    import pyretechnics.conversion as conv
    import pyretechnics.vector_utils as vu
    import pyretechnics.surface_fire as sf

@cy.cfunc
@cy.exceptval(check=False)
def van_wagner_critical_fireline_intensity(canopy_base_height: cy.float, foliar_moisture: cy.float) -> cy.float:
    """
    Returns the critical fireline intensity (kW/m) given:
    - canopy_base_height :: m
    - foliar_moisture    :: kg moisture/kg ovendry weight

    Constants used:
    460.0 = heat-of-ignition :: kJ/kg
    0.01 = empirical estimate for C in Van Wagner 1977 (eq. 4)
    """
    H: cy.float = 460.0 + 2600.0 * foliar_moisture
    v: cy.float = 0.01 * canopy_base_height * H
    return v * sqrt(v) # NOTE: This is faster than pow(v, 1.5).

@cy.cfunc
@cy.exceptval(check=False)
def van_wagner_crowning_spread_rate(surface_fire_max  : FireBehaviorMax,
                                    canopy_base_height: cy.float,
                                    foliar_moisture   : cy.float) -> cy.float:
    """
    Returns the surface spread rate above which crown fire occurs (m/min) given:
    - surface_fire_max   :: FireBehaviorMax struct
    - canopy_base_height :: m
    - foliar_moisture    :: kg moisture/kg ovendry weight
    """
    surface_max_fireline_intensity: cy.float = surface_fire_max.max_fireline_intensity
    if surface_max_fireline_intensity > 0.0:
        surface_max_spread_rate    : cy.float = surface_fire_max.max_spread_rate
        critical_fireline_intensity: cy.float = van_wagner_critical_fireline_intensity(canopy_base_height,
                                                                                       foliar_moisture)
        return (surface_max_spread_rate * critical_fireline_intensity / surface_max_fireline_intensity)
    else:
        return 0.0

@cy.ccall
@cy.inline
@cy.exceptval(check=False)
def van_wagner_crown_fire_initiation(surface_fireline_intensity: cy.float,
                                     canopy_cover              : cy.float,
                                     canopy_base_height        : cy.float,
                                     foliar_moisture           : cy.float) -> cy.bint:
    """
    Returns True if the surface fire transitions to a crown fire or False otherwise given:
    - surface_fireline_intensity :: kW/m
    - canopy_cover               :: 0-1
    - canopy_base_height         :: m
    - foliar_moisture            :: kg moisture/kg ovendry weight
    """
    return (
        surface_fireline_intensity > 0.0
        and
        canopy_cover > 0.4
        and
        surface_fireline_intensity >= van_wagner_critical_fireline_intensity(canopy_base_height, foliar_moisture)
    )

@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def cruz_active_crown_fire_spread_rate(wind_speed_10m              : cy.float,
                                       canopy_bulk_density         : cy.float,
                                       estimated_fine_fuel_moisture: cy.float) -> cy.float:
    """
    Returns the active crown fire spread rate (m/min) given:
    - wind_speed_10m                                   :: km/hr
    - canopy_bulk_density                              :: kg/m^3
    - estimated_fine_fuel_moisture (M_f[0] "dead-1hr") :: kg moisture/kg ovendry weight
    """
    return (11.02
            * pow(wind_speed_10m, 0.90)
            * pow(canopy_bulk_density, 0.19)
            * exp(-17.0 * estimated_fine_fuel_moisture))

@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def van_wagner_critical_spread_rate(canopy_bulk_density: cy.float) -> cy.float:
    """
    Returns the critical spread rate (m/min) given:
    - canopy_bulk_density :: kg/m^3
    """
    return 3.0 / canopy_bulk_density

@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def cruz_passive_crown_fire_spread_rate(active_spread_rate: cy.float, critical_spread_rate: cy.float) -> cy.float:
    """
    Returns the passive crown fire spread rate (m/min) given:
    - active_spread_rate   :: m/min
    - critical_spread_rate :: m/min
    """
    return active_spread_rate * exp(-active_spread_rate / critical_spread_rate)

@cy.cfunc
@cy.exceptval(check=False)
def cruz_crown_fire_spread_info(wind_speed_10m              : cy.float,
                                canopy_bulk_density         : cy.float,
                                estimated_fine_fuel_moisture: cy.float) -> CrownSpreadInfo:
    """
    Given these inputs:
    - wind_speed_10m                                   :: km/hr
    - canopy_bulk_density                              :: kg/m^3
    - estimated_fine_fuel_moisture (M_f[0] "dead-1hr") :: kg moisture/kg ovendry weight

    return a CrownSpreadInfo struct containing these keys:
    - fire_type            :: 2 (passive_crown) or 3 (active_crown)
    - spread_rate          :: m/min
    - critical_spread_rate :: m/min
    """
    active_spread_rate  : cy.float = cruz_active_crown_fire_spread_rate(wind_speed_10m,
                                                                        canopy_bulk_density,
                                                                        estimated_fine_fuel_moisture) # m/min
    critical_spread_rate: cy.float = van_wagner_critical_spread_rate(canopy_bulk_density) # m/min
    if (active_spread_rate > critical_spread_rate):
        return CrownSpreadInfo(
            fire_type            = 3, # FIXME: NAMED CONSTANT
            spread_rate          = active_spread_rate,
            critical_spread_rate = critical_spread_rate,
        )
    else:
        return CrownSpreadInfo(
            fire_type            = 2, # FIXME: NAMED CONSTANT
            spread_rate          = cruz_passive_crown_fire_spread_rate(active_spread_rate, critical_spread_rate),
            critical_spread_rate = critical_spread_rate,
        )

# NOTE: heat_of_combustion is h from the fuel models (generally 8000 Btu/lb).
# NOTE: ELMFIRE hard-codes heat_of_combustion to 18000 kJ/kg = 7738.6 Btu/lb.
@cy.cfunc
@cy.exceptval(check=False)
def calc_crown_fireline_intensity(crown_spread_rate  : cy.float,
                                  canopy_bulk_density: cy.float,
                                  canopy_height      : cy.float,
                                  canopy_base_height : cy.float,
                                  heat_of_combustion : cy.float) -> cy.float:
    """
    Returns the crown fireline intensity (Btu/ft/s OR kW/m) given:
    - crown_spread_rate   :: ft/min  OR m/min
    - canopy_bulk_density :: lb/ft^3 OR kg/m^3
    - canopy_height       :: ft      OR m
    - canopy_base_height  :: ft      OR m
    - heat_of_combustion  :: Btu/lb  OR kJ/kg

    (ft/min * lb/ft^3 * ft * Btu/lb)/60 = (Btu/ft/min)/60 = Btu/ft/s
    OR
    (m/min * kg/m^3 * m * kJ/kg)/60 = (kJ/m*min)/60 = kJ/m*s = kW/m
    """
    canopy_height_difference: cy.float = canopy_height - canopy_base_height
    return (crown_spread_rate * canopy_bulk_density * canopy_height_difference * heat_of_combustion) / 60.0

# Parameters for the linear model that computes LoW from wind speed.
@cy.cfunc
@cy.exceptval(check=False)
def crown_length_to_width_ratio(wind_speed_10m: cy.float, max_length_to_width_ratio: cy.float = 1e10) -> cy.float:
    """
    Calculate the length_to_width_ratio of the crown fire front using eq. 9 from
    Rothermel 1991 given:
    - wind_speed_10m            :: km/hr (aligned with the slope-tangential plane)
    - max_length_to_width_ratio :: float > 0.0 (Optional)
    """
    wind_speed_20ft      : cy.float = conv.wind_speed_10m_to_wind_speed_20ft(wind_speed_10m) # km/hr
    wind_speed_20ft_mph  : cy.float = conv.km_hr_to_mph(wind_speed_20ft)                     # mph
    length_to_width_ratio: cy.float = 1.0 + 0.125 * wind_speed_20ft_mph
    return min(length_to_width_ratio, max_length_to_width_ratio)


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def crown_fire_eccentricity(length_to_width_ratio: cy.float) -> cy.float:
    """
    Calculate the eccentricity (E) of the crown fire front using eq. 8 from
    Albini and Chase 1980 given:
    - length_to_width_ratio :: (1: circular spread, > 1: elliptical spread)
    """
    return sqrt(length_to_width_ratio * length_to_width_ratio - 1.0) / length_to_width_ratio

@cy.ccall
@cy.exceptval(check=False)
def calc_crown_fire_behavior_max(canopy_height               : cy.float,
                                 canopy_base_height          : cy.float,
                                 canopy_bulk_density         : cy.float,
                                 heat_of_combustion          : cy.float,
                                 estimated_fine_fuel_moisture: cy.float,
                                 wind_speed_10m              : cy.float,
                                 upwind_direction            : cy.float,
                                 slope                       : cy.float,
                                 aspect                      : cy.float,
                                 crown_max_lw_ratio          : cy.float=1e10) -> FireBehaviorMax:
    """
    Given these inputs:
    - canopy_height                                    :: m
    - canopy_base_height                               :: m
    - canopy_bulk_density                              :: kg/m^3
    - heat_of_combustion                               :: kJ/kg
    - estimated_fine_fuel_moisture (M_f[0] "dead-1hr") :: kg moisture/kg ovendry weight
    - wind_speed_10m                                   :: km/hr
    - upwind_direction                                 :: degrees clockwise from North
    - slope                                            :: rise/run
    - aspect                                           :: degrees clockwise from North
    - crown_max_lw_ratio                               :: float > 0.0 (Optional)

    return a FireBehaviorMax struct containing these keys:
    - max_fire_type          :: 2 (passive_crown) or 3 (active_crown)
    - max_spread_rate        :: m/min
    - max_spread_direction   :: (x, y, z) unit vector
    - max_fireline_intensity :: kW/m
    - max_flame_length       :: m
    - length_to_width_ratio  :: unitless (1: circular spread, > 1: elliptical spread)
    - eccentricity           :: unitless (0: circular spread, > 0: elliptical spread)
    - critical_spread_rate   :: m/min
    """
    # Reverse the provided wind and slope directions
    downwind_direction: cy.float = conv.opposite_direction(upwind_direction)
    upslope_direction : cy.float = conv.opposite_direction(aspect)
    # Project wind and slope vectors onto the slope-tangential plane
    # FIXME: Let's just have these vectors as arguments to the function instead of re-computing them.
    vectors: ProjectedVectors = sf.project_wind_and_slope_vectors_3d(wind_speed_10m, downwind_direction,
                                                                     slope, upslope_direction)
    wind_vector_3d : vec_xyz = vectors.wind_vector_3d  # km/hr
    slope_vector_3d: vec_xyz = vectors.slope_vector_3d # rise/run
    # Determine the max spread direction
    # FIXME: REVIEW Should we make the max_spread_direction the combined wind and slope direction?
    wind_speed_10m_3d: cy.float = vu.vector_magnitude_3d(wind_vector_3d) # km/hr
    max_spread_direction: vec_xyz
    if wind_speed_10m_3d > 0.0:
        max_spread_direction = vu.as_unit_vector_3d(wind_vector_3d)  # unit vector in the 3D downwind direction
    elif slope > 0.0:
        max_spread_direction = vu.as_unit_vector_3d(slope_vector_3d) # unit vector in the 3D upslope direction
    else:
        max_spread_direction = (0.0, 1.0, 0.0)                       # default: North
    # Calculate the crown fire behavior in the max spread direction
    spread_info          : CrownSpreadInfo = cruz_crown_fire_spread_info(wind_speed_10m_3d,
                                                                         canopy_bulk_density,
                                                                         estimated_fine_fuel_moisture)
    spread_rate          : cy.float        = spread_info.spread_rate                           # m/min
    fireline_intensity   : cy.float        = calc_crown_fireline_intensity(spread_rate,
                                                                           canopy_bulk_density,
                                                                           canopy_height,
                                                                           canopy_base_height,
                                                                           heat_of_combustion) # kW/m
    length_to_width_ratio: cy.float        = crown_length_to_width_ratio(wind_speed_10m_3d,
                                                                         crown_max_lw_ratio)   # unitless
    eccentricity         : cy.float        = crown_fire_eccentricity(length_to_width_ratio)    # unitless
    return FireBehaviorMax(
        max_fire_type          = spread_info.fire_type,
        max_spread_rate        = spread_rate,
        max_spread_direction   = max_spread_direction, # unit vector
        max_fireline_intensity = fireline_intensity,
        max_flame_length       = 0.0, # NOTE: max_flame_length is not provided, as in the original unoptimized code.
        length_to_width_ratio  = length_to_width_ratio,
        eccentricity           = eccentricity,
        critical_spread_rate   = spread_info.critical_spread_rate,
    )

@cy.ccall
@cy.exceptval(check=False)
def calc_crown_fire_behavior_in_direction(crown_fire_max  : FireBehaviorMax,
                                          spread_direction: vec_xyz) -> SpreadBehavior:
    """
    Given these inputs:
    - crown_fire_max     :: a FireBehaviorMax struct of max crown fire behavior values
      - max_fire_type          :: 2 (passive_crown) or 3 (active_crown)
      - max_spread_rate        :: m/min
      - max_spread_direction   :: (x, y, z) unit vector
      - max_fireline_intensity :: kW/m
      - max_flame_length       :: m
      - length_to_width_ratio  :: unitless (1: circular spread, > 1: elliptical spread)
      - eccentricity           :: unitless (0: circular spread, > 0: elliptical spread)
      - critical_spread_rate   :: m/min
    - spread_direction   :: 3D unit vector on the slope-tangential plane

    return a SpreadBehavior struct containing these keys:
    - dphi_dt            :: phi/min
    - fire_type          :: 2 (passive_crown) or 3 (active_crown)
    - spread_rate        :: m/min
    - spread_direction   :: (x, y, z) unit vector
    - fireline_intensity :: kW/m
    - flame_length       :: m
    """
    # Unpack max crown fire behavior values
    max_fire_type         : cy.int   = crown_fire_max.max_fire_type
    max_spread_rate       : cy.float = crown_fire_max.max_spread_rate
    max_spread_direction  : vec_xyz  = crown_fire_max.max_spread_direction
    max_fireline_intensity: cy.float = crown_fire_max.max_fireline_intensity
    eccentricity          : cy.float = crown_fire_max.eccentricity
    critical_spread_rate  : cy.float = crown_fire_max.critical_spread_rate
    # Calculate cos(w), where w is the offset angle between these unit vectors on the slope-tangential plane
    cos_w: cy.float = vu.dot_3d(max_spread_direction, spread_direction)
    # Calculate adjustment due to the offset angle from the max spread direction
    adjustment: cy.float = (1.0 - eccentricity) / (1.0 - eccentricity * cos_w)
    # Adjust the spread rate (possibly switching from an active to passive crown fire)
    spread_rate: cy.float = max_spread_rate * adjustment
    if spread_rate > critical_spread_rate:
        # Max spread rate was active and directional spread rate remains active
        return SpreadBehavior(
            dphi_dt            = 0.0,
            fire_type          = 3, # active_crown
            spread_rate        = spread_rate,
            spread_direction   = spread_direction,
            fireline_intensity = max_fireline_intensity * adjustment,
            flame_length       = 0.0,
        )
    elif max_fire_type == 2: # passive_crown
        # Max spread rate was passive and directional spread rate remains passive
        return SpreadBehavior(
            dphi_dt            = 0.0,
            fire_type          = 2, # passive_crown
            spread_rate        = spread_rate,
            spread_direction   = spread_direction,
            fireline_intensity = max_fireline_intensity * adjustment,
            flame_length       = 0.0,
        )
    else:
        # Max spread rate was active and directional spread rate has become passive
        return SpreadBehavior(
            dphi_dt            = 0.0,
            fire_type          = 2, # passive_crown
            spread_rate        = cruz_passive_crown_fire_spread_rate(spread_rate, critical_spread_rate),
            spread_direction   = spread_direction,
            fireline_intensity = max_fireline_intensity * adjustment,
            flame_length       = 0.0,
        )

from pprint import pprint
import pyretechnics.crown_fire as cf

# Set tree canopy characteristics
canopy_cover_low   = 0.3 # 0-1
canopy_cover_high  = 0.8 # 0-1
canopy_base_height = 3.0 # m
foliar_moisture    = 1.0 # kg moisture/kg ovendry weight

# Specify some surface fireline intensity values to check
# These would normally be calculated by the functions in pyretechnics.surface_fire
surface_fireline_intensity_low  = 500.0  # kW/m
surface_fireline_intensity_high = 5000.0 # kW/m

#=================================================================================
# Check for crown fire with canopy_cover_low and surface_fireline_intensity_low
#=================================================================================

print("Low Canopy Cover and Low Surface Fireline Intensity\n" + "=" * 60)

if cf.van_wagner_crown_fire_initiation(surface_fireline_intensity_low,
                                       canopy_cover_low,
                                       canopy_base_height,
                                       foliar_moisture):
    # Crown Fire
    print("\nA crown fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_low,
        "canopy_cover"              : canopy_cover_low,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })
else:
    # Surface Fire
    print("\nA surface fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_low,
        "canopy_cover"              : canopy_cover_low,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })

#=================================================================================
# Check for crown fire with canopy_cover_low and surface_fireline_intensity_high
#=================================================================================

print("\nLow Canopy Cover and High Surface Fireline Intensity\n" + "=" * 60)

if cf.van_wagner_crown_fire_initiation(surface_fireline_intensity_high,
                                       canopy_cover_low,
                                       canopy_base_height,
                                       foliar_moisture):
    # Crown Fire
    print("\nA crown fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_high,
        "canopy_cover"              : canopy_cover_low,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })
else:
    # Surface Fire
    print("\nA surface fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_high,
        "canopy_cover"              : canopy_cover_low,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })

#=================================================================================
# Check for crown fire with canopy_cover_high and surface_fireline_intensity_low
#=================================================================================

print("\nHigh Canopy Cover and Low Surface Fireline Intensity\n" + "=" * 60)

if cf.van_wagner_crown_fire_initiation(surface_fireline_intensity_low,
                                       canopy_cover_high,
                                       canopy_base_height,
                                       foliar_moisture):
    # Crown Fire
    print("\nA crown fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_low,
        "canopy_cover"              : canopy_cover_high,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })
else:
    # Surface Fire
    print("\nA surface fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_low,
        "canopy_cover"              : canopy_cover_high,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })

#=================================================================================
# Check for crown fire with canopy_cover_high and surface_fireline_intensity_high
#=================================================================================

print("\nHigh Canopy Cover and High Surface Fireline Intensity\n" + "=" * 60)

if cf.van_wagner_crown_fire_initiation(surface_fireline_intensity_high,
                                       canopy_cover_high,
                                       canopy_base_height,
                                       foliar_moisture):
    # Crown Fire
    print("\nA crown fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_high,
        "canopy_cover"              : canopy_cover_high,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })
else:
    # Surface Fire
    print("\nA surface fire occurs with these input parameters:")
    pprint({
        "surface_fireline_intensity": surface_fireline_intensity_high,
        "canopy_cover"              : canopy_cover_high,
        "canopy_base_height"        : canopy_base_height,
        "foliar_moisture"           : foliar_moisture,
    })

from pprint import pprint
import pyretechnics.crown_fire as cf

# Set input parameters
canopy_height                = 30.0    # m
canopy_base_height           = 3.0     # m
canopy_bulk_density          = 0.3     # kg/m^3
heat_of_combustion           = 18608.0 # kJ/kg
estimated_fine_fuel_moisture = 0.05    # kg moisture/kg ovendry weight
wind_speed_10m               = 10.0    # km/hr
upwind_direction             = 180.0   # degrees clockwise from North
slope                        = 0.5     # rise/run
aspect_parallel_to_wind      = 180.0   # degrees clockwise from North
aspect_perpendicular_to_wind = 270.0   # degrees clockwise from North

#====================================================================================================
# Calculate crown fire behavior in the direction of maximum spread with slope and wind aligned
#====================================================================================================

crown_fire_max_aligned = cf.calc_crown_fire_behavior_max(canopy_height, canopy_base_height,
                                                         canopy_bulk_density, heat_of_combustion,
                                                         estimated_fine_fuel_moisture,
                                                         wind_speed_10m, upwind_direction,
                                                         slope, aspect_parallel_to_wind)

# View results
print("Crown Fire Behavior Max with Wind and Slope Aligned\n" + "=" * 60)

print("\nInputs:")
pprint({
    "canopy_height"               : canopy_height,
    "canopy_base_height"          : canopy_base_height,
    "canopy_bulk_density"         : canopy_bulk_density,
    "heat_of_combustion"          : heat_of_combustion,
    "estimated_fine_fuel_moisture": estimated_fine_fuel_moisture,
    "wind_speed_10m"              : wind_speed_10m,
    "upwind_direction"            : upwind_direction,
    "slope"                       : slope,
    "aspect"                      : aspect_parallel_to_wind,
})

print("\nOutputs:")
pprint(crown_fire_max_aligned)

#====================================================================================================
# Calculate crown fire behavior in the direction of maximum spread with slope and wind perpendicular
#====================================================================================================

crown_fire_max_perpendicular = cf.calc_crown_fire_behavior_max(canopy_height, canopy_base_height,
                                                               canopy_bulk_density, heat_of_combustion,
                                                               estimated_fine_fuel_moisture,
                                                               wind_speed_10m, upwind_direction,
                                                               slope, aspect_perpendicular_to_wind)

# View results
print("\nCrown Fire Behavior Max with Wind and Slope Perpendicular\n" + "=" * 60)

print("\nInputs:")
pprint({
    "canopy_height"               : canopy_height,
    "canopy_base_height"          : canopy_base_height,
    "canopy_bulk_density"         : canopy_bulk_density,
    "heat_of_combustion"          : heat_of_combustion,
    "estimated_fine_fuel_moisture": estimated_fine_fuel_moisture,
    "wind_speed_10m"              : wind_speed_10m,
    "upwind_direction"            : upwind_direction,
    "slope"                       : slope,
    "aspect"                      : aspect_perpendicular_to_wind,
})

print("\nOutputs:")
pprint(crown_fire_max_perpendicular)

import numpy as np
from pprint import pprint
import pyretechnics.crown_fire as cf
import pyretechnics.vector_utils as vu

# Set input parameters
canopy_height                = 30.0    # m
canopy_base_height           = 3.0     # m
canopy_bulk_density          = 0.3     # kg/m^3
heat_of_combustion           = 18608.0 # kJ/kg
estimated_fine_fuel_moisture = 0.05    # kg moisture/kg ovendry weight
wind_speed_10m               = 10.0    # km/hr
upwind_direction             = 180.0   # degrees clockwise from North
slope                        = 0.5     # rise/run
aspect                       = 225.0   # degrees clockwise from North

#================================================================================================
# Calculate crown fire behavior in the direction of maximum spread
#================================================================================================

crown_fire_max = cf.calc_crown_fire_behavior_max(canopy_height, canopy_base_height,
                                                 canopy_bulk_density, heat_of_combustion,
                                                 estimated_fine_fuel_moisture,
                                                 wind_speed_10m, upwind_direction,
                                                 slope, aspect)

#================================================================================================
# Calculate crown fire behavior for the heading, flanking, and backing fires
#================================================================================================

heading_fire_direction  = crown_fire_max["max_spread_direction"]
# FIXME: Change upstream functions to return max_spread_direction as type vec_xyz
#        and change downstream functions to accept spread_direction as type vec_xyz
(dx, dy, dz)            = heading_fire_direction
flanking_fire_direction = np.asarray(vu.rotate_on_sloped_plane((dx, dy, dz), 90.0, slope, aspect))
backing_fire_direction  = np.asarray(vu.rotate_on_sloped_plane((dx, dy, dz), 180.0, slope, aspect))

heading_fire_behavior  = cf.calc_crown_fire_behavior_in_direction(crown_fire_max, heading_fire_direction)
flanking_fire_behavior = cf.calc_crown_fire_behavior_in_direction(crown_fire_max, flanking_fire_direction)
backing_fire_behavior  = cf.calc_crown_fire_behavior_in_direction(crown_fire_max, backing_fire_direction)

# View heading fire results
print("Heading Fire Behavior\n" + "=" * 60)
pprint(heading_fire_behavior)

# View flanking fire results
print("\nFlanking Fire Behavior\n" + "=" * 60)
pprint(flanking_fire_behavior)

# View backing fire results
print("\nBacking Fire Behavior\n" + "=" * 60)
pprint(backing_fire_behavior)

from pprint import pprint
import pyretechnics.conversion as conv
import pyretechnics.crown_fire as cf
import pyretechnics.fuel_models as fm
import pyretechnics.surface_fire as sf
import pyretechnics.vector_utils as vu

#================================================================================================
# Set input parameters
#================================================================================================

# Landscape
fuel_model          = fm.get_fuel_model(101) # GR1
canopy_cover        = 0.6                    # 0-1
canopy_height       = 30.0                   # m
canopy_base_height  = 3.0                    # m
canopy_bulk_density = 0.3                    # kg/m^3
slope               = 0.8                    # rise/run
aspect              = 225.0                  # degrees clockwise from North

# Weather
fuel_moisture    = (
    0.05, # dead_1hr
    0.10, # dead_10hr
    0.15, # dead_100hr
    0.00, # dead_herbaceous (will be set by moisturize for dynamic fuel models)
    0.90, # live_herbaceous
    0.60, # live_woody
)                        # kg moisture/kg ovendry weight
foliar_moisture  = 0.9   # kg moisture/kg ovendry weight
wind_speed_10m   = 30.0  # km/hr
upwind_direction = 180.0 # degrees clockwise from North

# Derived Parameters
fuel_bed_depth               = fuel_model["delta"]                      # ft
heat_of_combustion           = conv.Btu_lb_to_kJ_kg(fuel_model["h"][0]) # kJ/kg
estimated_fine_fuel_moisture = fuel_moisture[0]                         # kg moisture/kg ovendry weight

# Fire Perimeter Normal Vector
fire_perimeter_normal_vector_2d = (1.0, 0.0) # (x,y) unit vector pointing East

#================================================================================================
# Project the 2D fire perimeter normal vector onto the slope-tangential plane
#================================================================================================

upslope_direction               = conv.opposite_direction(aspect)
slope_vector_2d                 = conv.azimuthal_to_cartesian(slope, upslope_direction)
fire_perimeter_normal_vector_3d = vu.as_unit_vector_3d(vu.to_slope_plane(fire_perimeter_normal_vector_2d,
                                                                         slope_vector_2d))

#================================================================================================
# Calculate midflame wind speed
#================================================================================================

# Convert from 10m wind speed to 20ft wind speed
wind_speed_20ft = conv.wind_speed_10m_to_wind_speed_20ft(wind_speed_10m) # km/hr

# Convert 20ft wind speed from km/hr to m/min
wind_speed_20ft_m_min = conv.km_hr_to_m_min(wind_speed_20ft) # m/min

# Convert from 20ft wind speed to midflame wind speed in m/min
midflame_wind_speed = sf.calc_midflame_wind_speed(wind_speed_20ft_m_min,       # m/min
                                                  fuel_bed_depth,              # ft
                                                  conv.m_to_ft(canopy_height), # ft
                                                  canopy_cover)                # 0-1

#================================================================================================
# Calculate surface fire behavior in the direction of maximum spread
#================================================================================================

# Apply fuel moisture to fuel model
moisturized_fuel_model = fm.moisturize(fuel_model, fuel_moisture)

# Calculate no-wind-no-slope surface fire behavior
surface_fire_min = sf.calc_surface_fire_behavior_no_wind_no_slope(moisturized_fuel_model)

# Calculate surface fire behavior in the direction of maximum spread with effective wind speed limit
surface_fire_max = sf.calc_surface_fire_behavior_max(surface_fire_min,
                                                     midflame_wind_speed,
                                                     upwind_direction,
                                                     slope,
                                                     aspect,
                                                     surface_lw_ratio_model="rothermel")

#================================================================================================
# Calculate surface fire behavior in the direction normal to the fire perimeter
#================================================================================================

surface_fire_normal = sf.calc_surface_fire_behavior_in_direction(surface_fire_max, fire_perimeter_normal_vector_3d)

#================================================================================================
# Determine whether the surface fire transitions to a crown fire
#================================================================================================

if cf.van_wagner_crown_fire_initiation(surface_fire_normal["fireline_intensity"],
                                       canopy_cover,
                                       canopy_base_height,
                                       foliar_moisture):

    #============================================================================================
    # Calculate crown fire behavior in the direction of maximum spread
    #============================================================================================

    crown_fire_max = cf.calc_crown_fire_behavior_max(canopy_height, canopy_base_height,
                                                     canopy_bulk_density, heat_of_combustion,
                                                     estimated_fine_fuel_moisture,
                                                     wind_speed_10m, upwind_direction,
                                                     slope, aspect)

    #============================================================================================
    # Calculate crown fire behavior in the direction normal to the fire perimeter
    #============================================================================================

    crown_fire_normal = cf.calc_crown_fire_behavior_in_direction(crown_fire_max, fire_perimeter_normal_vector_3d)

    #============================================================================================
    # Calculate combined fire behavior in the direction normal to the fire perimeter
    #============================================================================================

    combined_fire_normal = cf.calc_combined_fire_behavior(surface_fire_normal, crown_fire_normal)

    #============================================================================================
    # Display the combined fire behavior normal to the fire perimeter
    #============================================================================================

    print("Fire Behavior from Combined Surface and Crown Fire\n" + "=" * 60)
    pprint(combined_fire_normal)

else:

    #============================================================================================
    # Display the surface fire behavior normal to the fire perimeter
    #============================================================================================

    print("Fire Behavior from Surface Fire (No Crown Fire Occurred)\n" + "=" * 60)
    pprint(surface_fire_normal)

import cython
import cython as cy
if cython.compiled:
    from cython.cimports.libc.math import round, sqrt, pow, log, exp, sin, cos
    from cython.cimports.pyretechnics.cy_types import \
        pyidx, vec_xy, coord_yx, coord_tyx, SpreadBehavior, SpotConfig, JumpDistribution
    from cython.cimports.pyretechnics.random import BufferedRandGen
    from cython.cimports.pyretechnics.space_time_cube import ISpaceTimeCube
    import cython.cimports.pyretechnics.conversion as conv
    import cython.cimports.pyretechnics.fuel_models as fm
    import cython.cimports.pyretechnics.surface_fire as sf
else:
    from math import sqrt, pow, log, exp, sin, cos
    from pyretechnics.py_types import \
        pyidx, vec_xy, coord_yx, coord_tyx, SpreadBehavior, SpotConfig, JumpDistribution
    from pyretechnics.random import BufferedRandGen
    from pyretechnics.space_time_cube import ISpaceTimeCube
    import pyretechnics.conversion as conv
    import pyretechnics.fuel_models as fm
    import pyretechnics.surface_fire as sf

@cy.cfunc
@cy.exceptval(check=False)
def expected_firebrand_production(fire_behavior           : SpreadBehavior,
                                  elevation_gradient      : vec_xy,
                                  cell_horizontal_area    : cy.float,
                                  firebrands_per_unit_heat: cy.float=1e-6) -> cy.float:
    """
    Return the expected number of firebrands produced by an entire cell when it burns given:
    - fire_behavior            :: a SpreadBehavior struct of surface or crown fire behavior values
      - dphi_dt                   :: phi/min
      - fire_type                 :: 0 (unburned), 1 (surface), 2 (passive_crown), or 3 (active_crown)
      - spread_rate               :: m/min
      - spread_direction          :: (x, y, z) unit vector on the slope-tangential plane
      - fireline_intensity        :: kW/m
      - flame_length              :: m
    - elevation_gradient       :: tuple with these fields
      - dz_dx                     :: rise/run
      - dz_dy                     :: rise/run
    - cell_horizontal_area     :: m^2
    - firebrands_per_unit_heat :: firebrands/kJ
    """
    if fire_behavior.spread_rate == 0.0:
        return 0.0
    else:
        #================================================================================================
        # Calculate the heat output per unit area
        #================================================================================================

        spread_rate         : cy.float = fire_behavior.spread_rate                                  # m/min
        fireline_intensity  : cy.float = fire_behavior.fireline_intensity                           # kW/m
        heat_output_per_area: cy.float = sf.calc_areal_heat_output(spread_rate, fireline_intensity) # kJ/m^2

        #================================================================================================
        # Calculate the slope-adjusted cell area
        #================================================================================================

        (dz_dx, dz_dy)         = elevation_gradient                            # (rise/run, rise/run)
        slope_factor: cy.float = sqrt(1.0 + (dz_dx * dz_dx) + (dz_dy * dz_dy)) # unitless
        cell_area   : cy.float = cell_horizontal_area * slope_factor           # m^2

        #================================================================================================
        # Calculate the expected number of firebrands produced in this cell
        #================================================================================================

        cell_heat_output: cy.float = heat_output_per_area * cell_area            # kJ
        firebrand_count : cy.float = firebrands_per_unit_heat * cell_heat_output # number of firebrands
        return firebrand_count

@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def delta_to_grid_dx(cos_wdir: cy.float, sin_wdir: cy.float, delta_x: cy.float, delta_y: cy.float) -> cy.float:
    """
    Computes the grid-aligned x coordinate of the delta vector, given the wind-aligned [ΔX ΔY] coordinates.
    Returns a signed distance (same unit as ΔX and ΔY).

    NOTE:
    - sin_wdir = wdir_x
    - cos_wdir = wdir_perp_x
    """
    return delta_x * sin_wdir + delta_y * cos_wdir


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def delta_to_grid_dy(cos_wdir: cy.float, sin_wdir: cy.float, delta_x: cy.float, delta_y: cy.float) -> cy.float:
    """
    Computes the grid-aligned y coordinate of the delta vector, given the wind-aligned [ΔX ΔY] coordinates.
    Returns a signed distance (same unit as ΔX and ΔY).

    NOTE:
    - cos_wdir = wdir_y
    - sin_wdir = wdir_perp_y
    """
    return delta_x * cos_wdir + delta_y * sin_wdir


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def distance_to_n_cells(distance: cy.float, cell_size: cy.float) -> cy.int:
    """
    Converts a delta expressed as a signed distance to one expressed as a number of grid cells.
    """
    return int(round(distance / cell_size))

@cy.cfunc
@cy.exceptval(check=False)
def resolve_exp_delta_x(spot_config: SpotConfig, fireline_intensity: cy.float, wind_speed_20ft: cy.float) -> cy.float:
    """
    Computes the expected value E[ΔX] (in meters) of the downwind spotting distance ΔX given:
    - spot_config        :: a SpotConfig struct of spotting parameters
      - random_seed                  :: seed for a numpy.random.Generator object
      - firebrands_per_unit_heat     :: firebrands/kJ
      - downwind_distance_mean       :: meters
      - fireline_intensity_exponent  :: downwind_distance_mean multiplier [I^fireline_intensity_exponent]
      - wind_speed_exponent          :: downwind_distance_mean multiplier [U^wind_speed_exponent]
      - downwind_variance_mean_ratio :: meters^2 / meter [downwind_variance_mean_ratio = Var(X) / E(X)]
      - crosswind_distance_stdev     :: meters
      - decay_distance               :: meters
    - fireline_intensity :: kW/m
    - wind_speed_20ft    :: m/s
    """
    downwind_distance_mean     : cy.float = spot_config.downwind_distance_mean
    fireline_intensity_exponent: cy.float = spot_config.fireline_intensity_exponent
    wind_speed_exponent        : cy.float = spot_config.wind_speed_exponent
    return (downwind_distance_mean
            * pow(fireline_intensity, fireline_intensity_exponent)
            * pow(wind_speed_20ft, wind_speed_exponent))


@cy.cfunc
@cy.exceptval(check=False)
def resolve_var_delta_x(spot_config: SpotConfig, exp_delta_x: cy.float) -> cy.float:
    """
    Computes the variance Var[ΔX] (in m^2) of the downwind spotting distance ΔX given:
    - spot_config :: a SpotConfig struct of spotting parameters
      - random_seed                  :: seed for a numpy.random.Generator object
      - firebrands_per_unit_heat     :: firebrands/kJ
      - downwind_distance_mean       :: meters
      - fireline_intensity_exponent  :: downwind_distance_mean multiplier [I^fireline_intensity_exponent]
      - wind_speed_exponent          :: downwind_distance_mean multiplier [U^wind_speed_exponent]
      - downwind_variance_mean_ratio :: meters^2 / meter [downwind_variance_mean_ratio = Var(X) / E(X)]
      - crosswind_distance_stdev     :: meters
      - decay_distance               :: meters
    - exp_delta_x :: meters (E[ΔX])
    """
    downwind_variance_mean_ratio: cy.float = spot_config.downwind_variance_mean_ratio
    return downwind_variance_mean_ratio * exp_delta_x


@cy.cfunc
@cy.exceptval(check=False)
def lognormal_mu_from_moments(mean: cy.float, variance: cy.float) -> cy.float:
    """
    TODO: Add docstring
    """
    m2: cy.float = mean * mean
    return log(m2 / sqrt(m2 + variance))


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def lognormal_sigma_from_moments(mean: cy.float, variance: cy.float) -> cy.float:
    """
    TODO: Add docstring
    """
    return sqrt(log(1.0 + variance / (mean * mean)))


@cy.cfunc
@cy.exceptval(check=False)
def resolve_lognormal_params(spot_config       : SpotConfig,
                             fireline_intensity: cy.float,
                             wind_speed_20ft   : cy.float) -> tuple[cy.float, cy.float]:
    """
    TODO: Add docstring
    """
    exp_delta_x         : cy.float = resolve_exp_delta_x(spot_config, fireline_intensity, wind_speed_20ft)
    var_delta_x         : cy.float = resolve_var_delta_x(spot_config, exp_delta_x)
    prob_lognormal_mu   : cy.float = lognormal_mu_from_moments(exp_delta_x, var_delta_x)
    prob_lognormal_sigma: cy.float = lognormal_sigma_from_moments(exp_delta_x, var_delta_x)
    return (prob_lognormal_mu, prob_lognormal_sigma)

# When will we have the default sigma_Y > E[ΔX]?
# It can be seen that this nonsensical situation
# happens iff sigma_X exceeds the following number:
#
# sqrt(log(1.0 + (0.88 ** 2.0) / (0.92 * 0.47))
#
# => 1.0131023746492023
sigma_y_scalar_m = cy.declare(cy.double, 0.92 * 0.47 / (0.88 * 0.88))


@cy.cfunc
@cy.exceptval(check=False)
def himoto_resolve_default_sigma_y_from_lognormal_params(mu_x: cy.float, sigma_x: cy.float) -> cy.float:
    es2h      : cy.float = exp((sigma_x * sigma_x) / 2.0)
    avg_deltax: cy.float = exp(mu_x) * es2h
    return sigma_y_scalar_m * avg_deltax * (es2h + 1.0) * (es2h - 1.0) # meters


@cy.cfunc
@cy.exceptval(check=False)
def himoto_resolve_default_sigma_y(spot_config       : SpotConfig,
                                   fireline_intensity: cy.float,
                                   wind_speed_20ft   : cy.float) -> cy.float:
    ln_params: tuple[cy.float, cy.float] = resolve_lognormal_params(spot_config, fireline_intensity, wind_speed_20ft)
    mu_x     : cy.float                  = ln_params[0] # dimensionless (log-space)
    sigma_x  : cy.float                  = ln_params[1] # dimensionless (log-space)
    return himoto_resolve_default_sigma_y_from_lognormal_params(mu_x, sigma_x) # meters


@cy.cfunc
@cy.exceptval(check=False)
def resolve_crosswind_distance_stdev(spot_config       : SpotConfig,
                                     fireline_intensity: cy.float,
                                     wind_speed_20ft   : cy.float) -> cy.float:
    crosswind_distance_stdev: cy.float = spot_config.crosswind_distance_stdev
    if crosswind_distance_stdev != 0.0:
        return crosswind_distance_stdev # meters
    else:
        return himoto_resolve_default_sigma_y(spot_config, fireline_intensity, wind_speed_20ft) # meters

@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def sample_normal(rng: BufferedRandGen, mu: cy.float, sd: cy.float) -> cy.float:
    """
    Returns sample from normal/gaussian distribution given mu and sd.
    """
    return mu + sd * rng.next_normal()


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def sample_lognormal(rng: BufferedRandGen, mu: cy.float, sd: cy.float) -> cy.float:
    """
    Returns sample from log-normal distribution given mu and sd.
    """
    return exp(sample_normal(rng, mu, sd))


@cy.cfunc
@cy.exceptval(check=False)
def resolve_JumpDistribution(spot_config       : SpotConfig,
                             fireline_intensity: cy.float,
                             wind_speed_20ft   : cy.float) -> JumpDistribution:
    ln_params: tuple[cy.float, cy.float] = resolve_lognormal_params(spot_config, fireline_intensity, wind_speed_20ft)
    # Initialize a new JumpDistribution
    return JumpDistribution(
        mu_x    = ln_params[0], # dimensionless (log-space)
        sigma_x = ln_params[1], # dimensionless (log-space)
        sigma_y = resolve_crosswind_distance_stdev(spot_config, fireline_intensity, wind_speed_20ft), # meters
    )


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def sample_downwind_jump(jd: JumpDistribution, random_generator: BufferedRandGen) -> cy.float:
    return sample_lognormal(random_generator, jd.mu_x, jd.sigma_x)


@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def sample_crosswind_jump(jd: JumpDistribution, random_generator: BufferedRandGen) -> cy.float:
    return sample_normal(random_generator, 0.0, jd.sigma_y)

@cy.cfunc
@cy.exceptval(check=False)
def heat_of_preignition(temperature: cy.float, fine_fuel_moisture: cy.float) -> cy.float:
    """
    Returns heat of preignition given:
    - temperature        :: degrees Celsius
    - fine_fuel_moisture :: 0-1

    Q_ig = 144.512 - 0.266 * T_o - 0.00058 * (T_o)^2 - T_o * M + 18.54 * (1 - exp(-15.1 * M)) + 640 * M (eq. 10)
    """
    T_o: cy.float = temperature
    M  : cy.float = fine_fuel_moisture
    # Heat required to reach ignition temperature
    Q_a: cy.float = 144.512 - 0.266 * T_o - 0.00058 * (T_o * T_o)
    # Heat required to raise moisture to the boiling point
    Q_b: cy.float = -T_o * M
    # Heat of desorption
    Q_c: cy.float = 18.54 * (1.0 - exp(-15.1 * M))
    # Heat required to vaporize moisture
    Q_d: cy.float = 640.0 * M
    return Q_a + Q_b + Q_c + Q_d


@cy.ccall
@cy.exceptval(check=False)
def schroeder_ignition_probability(temperature: cy.float, fine_fuel_moisture: cy.float) -> cy.float:
    """
    Returns the probability of spot fire ignition (Schroeder 1969) given:
    - temperature        :: degrees Celsius
    - fine_fuel_moisture :: 0-1

    X    = (400 - Q_ig) / 10
    P(I) = (0.000048 * X^4.3) / 50 (pg. 15)
    """
    Q_ig: cy.float = heat_of_preignition(temperature, fine_fuel_moisture)
    X   : cy.float = max(0.0, 400.0 - Q_ig) * 0.1
    P_I : cy.float = 0.000048 * pow(X, 4.3) * 0.02
    return min(P_I, 1.0)

@cy.ccall
@cy.inline
@cy.exceptval(check=False)
def firebrand_flight_survival_probability(spotting_distance: cy.float, decay_distance: cy.float) -> cy.float:
    """
    Returns the probability that a firebrand will survive its flight (Perryman 2012) given:
    - spotting_distance :: meters (d)
    - decay_distance    :: meters (1/lambda)

    P(Survival) = exp(-d * lambda)
    """
    return exp(-spotting_distance / decay_distance)

# FIXME: unused
@cy.cfunc
@cy.inline
@cy.exceptval(check=False)
def albini_firebrand_maximum_height(firebrand_diameter: cy.float) -> cy.float:
    return 0.39e5 * firebrand_diameter


@cy.cfunc
@cy.exceptval(check=False)
def albini_t_max(flame_length: cy.float) -> cy.float:
    """
    Returns the time of spot ignition using Albini1979spot in minutes given:
    - flame_length :: meters [z_F]

    a           = 5.963                                                             (D33)
    b           = a - 1.4                                                           (D34)
    D           = 0.003
    z           = 0.39 * D * 10^5
    w_F         = 2.3 * z_F^0.5                                                     (A58)
    t_c         = 1
    t_o         = t_c / (2 * z_F / w_F)
    travel_time = t_1 + t_2 + t_3 = 1.2 + (a / 3) * (((b + (z/z_F)) / a)^3/2 - 1)   (D43)
    """
    a        : cy.float = 5.963  # dimensionless constant from (D33)
    b        : cy.float = 4.563  # dimensionless constant from (D34)
    z        : cy.float = 117.0  # maximum altitude of firebrands in meters [derived for (D44) in (Albini1979spot)]
    z_F      : cy.float = flame_length                      # m
    w_F      : cy.float = 2.3 * sqrt(flame_length)          # m/s
    charact_t: cy.float = conv.sec_to_min(2.0 * z_F / w_F)  # min
    # The following dimensionless factor is equal to t_T - t_o, with t_T defined by (D43) in Albini1979spot.
    u          : cy.float = (b + z / z_F) / a
    u3_2       : cy.float = u * sqrt(u) # Faster than ** 1.5
    travel_time: cy.float = 1.2 + (a / 3.0) * (u3_2 - 1.0)
    return charact_t * travel_time


# FIXME: Consider removing t_max from this calculation for performance.
@cy.cfunc
@cy.exceptval(check=False)
def spot_ignition_time(time_of_arrival: cy.float, flame_length: cy.float) -> cy.float:
    """
    Returns the time of spot ignition using Albini 1979 and Perryman 2012 in minutes given:
    - time_of_arrival :: minutes
    - flame_length    :: meters

    t_spot = time_of_arrival + (2 * t_max) + t_ss
    """
    t_max         : cy.float = albini_t_max(flame_length)
    t_steady_state: cy.float = 20.0 # period of building up to steady state from ignition (min)
    return time_of_arrival + 2.0 * t_max + t_steady_state