diff --git a/src/crop/new-almass.jl b/src/crop/new-almass.jl new file mode 100644 index 0000000000000000000000000000000000000000..31c52f3e8ca544e86fa71fee73656245064d8704 --- /dev/null +++ b/src/crop/new-almass.jl @@ -0,0 +1,261 @@ +# 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