Skip to content
Snippets Groups Projects
Commit 9184287d authored by Marco Matthies's avatar Marco Matthies
Browse files

Add crop model functions directly translated from ALMaSS

parent 892abe23
No related branches found
No related tags found
No related merge requests found
# 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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment