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
# - almass_RecalculateBugsNStuff()
# - tov_MaizeStrigling, etc
# - radconv = model.SolarConversion[is_c4_plant + 1][index + 1]


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

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

    newgrowth::Float64 = 0
    veg_cover::Float64 = 0

    vege_type::Int = 1
    oldLAtotal::Float64 = 0
    biomass_scale::Vector{Float64} = [1.0]
    owner_index::Int = 1
    poly::Int = 1
    veg_biomass::Float64 = 0
    area::Float64 = 1
    total_biomass::Float64 = 0
    total_biomass_old::Float64 = 0
    veg_density::Float64 = 0
    green_biomass::Float64 = 0
    dead_biomass::Float64 = 0
end

function Base.show(io::IO, ve::VegElement)
    println(io, "VegElement:")
    for field in fieldnames(typeof(ve))
        println(io, "    $field: ", getfield(ve, field))
    end
end

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

# `class CropGrowth` from `Landscape/Plants.h`, line 67
@kwdef struct CropGrowth
    # Note:
    # - 5 growth phases (enum Growth_Phases)
    # - 3 curves (leaf_area_green = 0, leaf_area_total, height)
    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

        pgd = get_plant_growth_data(model)
        dLAG = almass_GetLAgreenDiffScaled(pgd, ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
        dLAT = almass_GetLAtotalDiffScaled(pgd, ve.ddegs, ve.yddegs, ve.curve_num, ve.veg_phase, ve.growth_scaler)
        dHgt = almass_GetHeightDiffScaled(pgd, 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

    # Note: this calculates the biomass etc, not just "BugsNStuff"
    almass_RecalculateBugsNStuff(ve, model)

    # Note: missing code
    # - # 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

# `double Weather::GetDDDegs(long a_date)` from `Landscape/Weather.cpp`, line 205
function almass_GetDDDegs(model::SimulationModel)
    basetemp = 0.0  # TODO: doesn't seem to use a min growth temp ???
    #return bounds((maxtemp(model) + mintemp(model)) / 2 - basetemp)
    return max((maxtemp(model) + mintemp(model)) / 2 - basetemp, 0.0)
end

# `double PlantGrowthData::GetLAgreenDiffScaled(...)` from `Landscape/Plants.h`, line 136
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

# `double PlantGrowthData::GetLAtotalDiffScaled(...)` from `Landscape/Plants.h`, line 138
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

# `double PlantGrowthData::GetHeightDiffScaled(...)` from `Landscape/Plants.h`, line 140
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

# `double PlantGrowthData::GetLAgreenDiff(...)` from `Landscape/Plants.h`, line 106
function almass_GetLAgreenDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64,
                               curve_num::Int, veg_phase::Int)
    # TODO: hardcoded 1, in original code it is 0
    return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 1)
end

# `double PlantGrowthData::GetLAtotalDiff(...)` from `Landscape/Plants.h`, line 116
function almass_GetLAtotalDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64,
                               curve_num::Int, veg_phase::Int)
    # TODO: hardcoded 2, in original code it is 1
    return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 2)
end

# `double PlantGrowthData::GetHeightDiff(...)` from `Landscape/Plants.h`, line 125
function almass_GetHeightDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64,
                              curve_num::Int, veg_phase::Int)
    # TODO: hardcoded 3, in original code it is 2
    return almass_FindDiff(pgd, ddegs, yddegs, curve_num, veg_phase, 3)
end

# `double PlantGrowthData::FindDiff(...)` from `Landscape/Plants.cpp`, line 45
function almass_FindDiff(pgd::PlantGrowthData, ddegs::Float64, yddegs::Float64,
                         curve_num::Int, veg_phase::Int, type::Int)
    # Note: curve_num is a_plant in original code
    #
    # 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

# `void VegElement::RecalculateBugsNStuff(void)` from `Landscape/Elements.cpp`, line 1704
function almass_RecalculateBugsNStuff(ve::VegElement, model::SimulationModel)
    # Calculate vegetation cover and vegetation biomass
    #
    # Original comment:
    #
    # /** This is the heart of the dynamics of vegetation elements. It
    # calculates vegetation cover and uses this to determine
    # vegetation biomass.  It also calculates spilled grain and goose
    # forage, as well a calculating insect biomass, vegetation density
    # and dead biomass*/

    ve.newgrowth = 0
    # `cfg_beer_law_extinction_coef` from `Landscape/Elements.cpp`, line 139
    beer_law_extinction_coef = 0.6
    # Beer's Law to give cover
    ve.veg_cover = 1.0 - exp(- beer_law_extinction_coef * ve.LAtotal)
    # This is used to calc growth rate
    usefull_veg_cover = 1.0 - exp(-0.4 * ve.LAgreen)

    # global radiation today
    glrad = almass_SupplyGlobalRadiation(model)

    # This is different for maize (a C4 plant)
    is_c4_plant = ((ve.vege_type == tov_Maize) ||
        (ve.vege_type == tov_OMaizeSilage) ||
        (ve.vege_type == tov_MaizeSilage) ||
        (ve.vege_type == tov_MaizeStrigling))
    # TODO: simplify calculation of `index`
    index = Int(floor(0.5 + almass_SupplyTemp(model)) + 30 + 1) # // There are 30 negative temps
    # TODO: model needs this field???
    # Note: +1 because Julia arrays are 1-based
    radconv = almass_SolarConversion(model, is_c4_plant, index)

    if ve.LAtotal >= ve.oldLAtotal
        # we are in positive growth so grow depending on our equation
        ve.newgrowth = usefull_veg_cover * glrad * radconv * ve.biomass_scale[ve.vege_type];
        if ve.owner_index != -1
            # This only works because only crops and similar structures have owners
            fintensity = almass_SupplyFarmIntensity(model, ve.poly)
            if fintensity >= 1
                # 1 means extensive, so reduce vegetation biomass by 20%
                # NB this cannot be used with extensive crop types otherwise you get an additional 20% reduction
                # This way of doing things provides a quick and dirty general effect.
                ve.veg_biomass += ve.newgrowth * 0.8
            else
                ve.veg_biomass += ve.newgrowth
            end
        else
            ve.veg_biomass += ve.newgrowth
        end
    else
        # Negative growth - so shrink proportional to the loss in LAI Total
        if ve.oldLAtotal > 0
            temp_propotion = ve.LAtotal / ve.oldLAtotal
            # //m_newgrowth = m_veg_biomass*(temp_propotion-1);
            ve.veg_biomass *= temp_propotion;
        end
    end

    # Here we also want to know how much biomass we have on the field
    # in total. So we multiply the current biomass by area
    ve.total_biomass = ve.veg_biomass * ve.area
    ve.total_biomass_old = ve.total_biomass
    # NB The m_weed_biomass is calculated directly from the curve in Curves.pre
    # rather than going through the rigmarole of converting leaf-area index
    ve.veg_density = Int(floor(0.5 + (ve.veg_biomass / (1 + ve.veg_height))))
    if ve.veg_density > 100
        # to stop array bounds problems
        ve.veg_density = 100
    end

    if ve.LAtotal == 0.0
        ve.green_biomass = 0.0
    else
        ve.green_biomass = ve.veg_biomass * (ve.LAgreen / ve.LAtotal)
    end
    ve.dead_biomass = ve.veg_biomass - ve.green_biomass
    #// Here we use our member function pointer to call the correct piece of code for our current species
    #(this->*(SpeciesSpecificCalculations))();
end