new-almass.jl

# Functions translated from ALMaSS for use in Persefone
#
# ALMaSS is licensed as follows:
#
# Copyright (c) 2011, Christopher John Topping, Aarhus University
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided
# that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list of conditions and the
# following disclaimer.
# Redistributions in binary form must reproduce the above copyright notice, this list of conditions and
# the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
# FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
# BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


# TODO

# `class VegElement : public LE` from `Landscape/Elements.h`, line 601
struct VegElement
    # TODO: this class has a lot of fields, only the fields used are
    # here

    # new_weed_growth::Float64
    force_growth::Bool
    yddegs::Float64
    ddegs::Float64
    vegddegs::Float64
    curve_num::Int
    veg_phase::Int  # TODO: enum type?
    growth_scaler::Float64
    LAgreen::Float64
    LAtotal::Float64
    veg_height::Float64
    force_LAgreen::Float64
    force_LAtotal::Float64
    force_veg_height::Float64
end

# `typedef enum { ... } Growth_Phases` from `Landscape/Plants.h`, line 51
@enum Growth_Phases begin
    janfirst = 0
    sow
    marchfirst
    harvest1
    harvest2
end

# `class CropGrowth` from `Landscape/Plants.h`, line 67
# Note:
# - 5 growth phases (enum Growth_Phases)
# - 3 curves (leaf_area_green = 0, leaf_area_total, height)
struct CropGrowth
    lownut::Bool  # = false
    dds::Vector{Vector{Float64}}  # ddegs at inflection, double m_dds[5][MaxNoInflections] = { 0 }
    slopes::Vector{Vector{Vector{Float64}}}  # double m_slopes[5][3][MaxNoInflections] = { 0 }
    start::Vector{Vector{Float64}}  # double m_start[5][3] = { 0 }
    start_valid::Vector{Bool}  # bool m_start_valid[5] = { 0 }
end

# `class PlantGrowthData` from `Landscape/Plants.h`, line 79
struct PlantGrowthData
    growth::Vector{CropGrowth}
    final_ddeg::Vector{Vector{Vector{Float64}}}
    numbers::Vector{Int}
    num_crops::Int
    # Note: some more fields ommitted related to weeds and bugs
end

# `void VegElement::DoDevelopment(void)` from `Landscape/Elements.cpp`, line 2254
function almass_DoDevelopment(ve::VegElement, model::SimulationModel)
    # ve.new_weed_growth = 0.0
    if ! ve.force.growth
        # first do the day degree calculations
        ve.yddegs = ve.ddegs
        pos_temp_today = almass_GetDDDegs(model)
        if ve.vegddegs != -1.0
            # sum up the vegetation day degrees since sowing
            ve.vegddegs += pos_temp_today
        end
        # sum up the phase day degrees
        ve.ddegs = pos_temp_today + ve.yddegs

        dLAG = almass_GetLAgreenDiffScaled(ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
        dLAT = almass_GetLAtotalDiffScaled(ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
        dHgt = almass_GetHeightDiffScaled(ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)

        ve.LAgreen += dLAG
        if ve.LAgreen < 0.0
            ve.LAgreen = 0.0
        end
        ve.LAtotal += dLAT
        if ve.LAtotal < 0.0
            ve.LAtotal = 0.0
        end
        ve.veg_height += dHgt
        if ve.veg_height < 0.0
            ve.veg_height = 0.0
        end

        # TODO: ignored code for weeds
        #       if (this->m_owner_index != -1) ...
        #           // This only works because only crops and similar structures have owners
    else
        almass_ForceGrowthDevelopment(ve)
    end

    # check that ve.LAtotal is bigger than ve.LAgreen
    if ve.LAtotal < ve.LAgreen
        @warn "ve.LAtotal < ve.LAgreen, $(ve.LAtotal) < $(ve.LAgreen). Performing hack correction"
        ve.LAtotal = 1.1 * ve.LAgreen
    end

    # TODO: missing code
    # - RecalculateBugsNStuff()
    # - # Here we need to set today's goose numbers to zero in case they are not written
    #   # by the goose population manager (the normal situation)
    #   ResetGeese()
    # - # Deal with any possible unsprayed margin, transferring info as necessary
    #   if (GetUnsprayedMarginPolyRef() != -1) ...
    # - m_interested_green_biomass = m_green_biomass * m_interested_biomass_fraction
end

# TODO: original source location
function almass_GetDDDegs(model::SimulationModel)
    # Original code
    #
    # double temp = m_temp[ a_date%m_datemodulus ];
    # if ( temp < 0.0 ) {
    #     temp = 0.0;
    # }
    # return temp;

    basetemp = 0.0  # TODO: doesn't seem to use a min growth temp ???
    return bounds((maxtemp(model) + mintemp(model)) / 2 - basetemp)
end


# TODO: original source location
function almass_GetLAgreenDiffScaled(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64, curve_num::Int, veg_phase::Int, growth_scaler::Float64)
    return growth_scaler * almass_GetLAgreenDiff(pgd, ddegs, yddegs, curve_num, veg_phase)
end

# TODO: original source location
function almass_GetLAtotalDiffScaled(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64, curve_num::Int, veg_phase::Int, growth_scaler::Float64)
    return growth_scaler * almass_GetLAtotalDiff(pgd, ddegs, yddegs, curve_num, veg_phase)
end

# TODO: original source location
function almass_GetHeightDiffScaled(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64, curve_num::Int, veg_phase::Int, growth_scaler::Float64)
    return growth_scaler * almass_GetHeightDiff(pgd, ddegs, yddegs, curve_num, veg_phase)
end

# TODO: original source location
function almass_GetLAgreenDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64, curve_num::Int, veg_phase::Int)
    # TODO: hardcoded 0
    return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 0)
end

# TODO: original source location
function almass_GetLAtotalDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64, curve_num::Int, veg_phase::Int)
    # TODO: hardcoded 1
    return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 1)
end

function almass_GetHeightDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64, curve_num::Int, veg_phase::Int)
    # TODO: hardcoded 2
    return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 2)
end

# `double PlantGrowthData::FindDiff(...)` from `Landscape/Plants.cpp`, line 45
# Note: curve_num is a_plant in original code
function almass_FindDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64, curve_num::Int, veg_phase::Int, type::Int)
    # Note: curve_num is an index like in the CropGrowth.txt file,
    #       *not* indices into the m_growth array

    # TODO: Comments from original code
    #
    # // Check for valid plant number at runtime?
    # // This is broken for growth curves where one can risk passing
    # // more than a single inflection point in the growth curve in a
    # // single day...

    index = pgd.numbers[curve_num]  # index = m_numbers[ a_plant ]
    oldindex = 0
    newindex = 0

    # TODO: magic constant 99999
    if first(pgd.growth[index].dds[veg_phase]) == 99999
        return 0.0
    end

    dds_for_phase = pgd.growth[index].dds[veg_phase]

    # TODO: breaks on last index
    for i in eachindex(dds_for_phase)
        # // In other words: If the current value for summed day degrees
        # // is smaller than the X position of the *next* inflection
        # // point, then we are in the correct interval.
        if i == lastindex(dds_for_phase) || dds_for_phase[i + 1] > ddegs
            newindex = i
            break
        end
    end

    # TODO: illegal access for last index
    for i in eachindex(dds_for_phase)
        if i == lastindex(dds_for_phase) || dds_for_phase[i + 1] > yddegs
            oldindex = i
            break
        end
    end

    if newindex > oldindex
        # // We have passed an inflection point between today and yesterday.
        # // First add the increment from yesterdays day degree sum up to
        # // the inflection point.
        dddif = dds_for_phase[newindex] - yddegs
        diff = pgd.growth[index].slopes[veg_phase][type][oldindex] * dddif

        # // Then from the inflection point up to today.
        dddif = ddegs - dds_for_phase[newindex]
        diff += pgd.growth[index].slopes[veg_phase][type][newindex] * dddif
    else
        # // No inflection point passed.
        diff = pgd.growth[index].slopes[veg_phase][type][newindex] * (ddegs - yddegs)
    end

    return diff
end

# `double VegElement::ForceGrowthDevelopment(...)` from `Landscape/Elements.cpp`, line 2223
function almass_ForceGrowthDevelopment(ve::VegElement)
    # if ( m_herbicidedelay <= 0 ) {
    #     m_weed_biomass += m_force_Weed; // ***CJT*** 12th Sept 2008 - rather than force growth, weeds might be allowed to grow on their own
    #   if(m_force_Weed>0) m_new_weed_growth = m_force_Weed;
    # }

    ve.LAgreen += ve.force_LAgreen
    ve.LAtotal += ve.force_LAtotal
    ve.veg_height += ve.force_veg_height

    if ve.LAgreen < 0
        ve.LAgreen = 0
    end
    if ve.LAtotal < 0
        ve.LAtotal = 0
    end
    if ve.veg_height < 0
        ve.veg_height = 0
    end
end