plants.cpp

Go to the documentation of this file.

//
// plants.cpp
//
/*
*******************************************************************************************************
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.
********************************************************************************************************
*/
 
#define _CRT_SECURE_NO_DEPRECATE
 
#include <cstdio>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include "configurator.h"
#include "maperrormsg.h"
#include "ls.h"
 
extern CfgBool cfg_pollen_nectar_on;
extern void FloatToDouble(double&, float);
class PlantGrowthData *g_crops;
 
double PlantGrowthData::FindDiff( double a_ddegs, double a_yddegs, int a_plant,
                          int a_phase, int a_type )
{
  // 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...
 
  int index = m_numbers[ a_plant ];
  unsigned int oldindex=0, newindex=0;
 
  if ( m_growth[ index ]->m_dds[ a_phase ][ 0 ] == 99999 ) {
    return 0.0;
  }
 
  for ( unsigned int i=0; i<MaxNoInflections; i++ ) {
    // 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 ( m_growth[ index ]->m_dds[ a_phase ][ i+1 ] > a_ddegs ) {
      newindex = i;
      break;
      //      return m_growth[ index ]->m_slopes[ a_phase ][ a_type ][i];
    }
  }
 
  for ( unsigned int i=0; i<MaxNoInflections; i++ ) {
    if ( m_growth[ index ]->m_dds[ a_phase ][ i+1 ] > a_yddegs ) {
      oldindex = i;
      break;
      //      return m_growth[ index ]->m_slopes[ a_phase ][ a_type ][i];
    }
  }
 
  double diff;
 
  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.
    double dddif =
      m_growth[ index ]->m_dds[ a_phase ][ newindex ] - a_yddegs;
    diff =
      m_growth[ index ]->m_slopes[ a_phase ][ a_type ][oldindex]*
      dddif;
 
    // Then from the inflection point up to today.
    dddif = a_ddegs -
      m_growth[ index ]->m_dds[ a_phase ][ newindex ];
    diff +=
      m_growth[ index ]->m_slopes[ a_phase ][ a_type ][ newindex ]*
      dddif;
  } else {
    // No inflection point passed.
    diff = m_growth[ index ]->m_slopes[ a_phase ][ a_type ][ newindex ] *
      (a_ddegs - a_yddegs);
  }
  return diff;
}
 
unsigned int PlantGrowthData::FindCropNum( ifstream& ist )
{
  int NoPlants;
 
  m_numbers.resize(201);
  for ( unsigned int i=0; i<201; i++) {
    m_numbers[ i ] = -1;
  }
 
  //m_ifile = fopen(a_cropcurvefile, "r" );
  //ifstream m_ifile(a_cropcurvefile);
  /*
  if (!m_ifile){
    g_msg->Warn(WARN_FILE, "PlantGrowthData::PlantGrowthData: Unable to open file",
                a_cropcurvefile );
    exit(1);
  }
  */
  //fscanf( m_ifile , "%d", &NoPlants );  // How many tables to read in
  ist >> NoPlants;
  m_growth.resize( NoPlants );
  m_num_crops = NoPlants;
 
  return NoPlants;
}
 
 
 
void PlantGrowthData::SetVegNum( unsigned int a_i, ifstream& ist, const char * a_cropcurvefile )
{
  int ThisPlant;
 
  // Find out what crop and what nutrient status
  //fscanf( m_ifile, "%d", &ThisPlant);
 
  ist >> ThisPlant;
  // Check if valid plant number (from the file).
  if ( ThisPlant < 0 || ThisPlant > 200 ) {
    g_msg->Warn(WARN_FILE, "PlantGrowthData::FindCropNum(): Illegal plant number"
                " specified in", a_cropcurvefile );
    exit(1);
  }
 
  m_numbers[ ThisPlant ] = a_i;
 
  // if greater than 100 then it is low nutrient
  if ( ThisPlant > 100 ) {
    m_growth[ a_i ]->m_lownut = true;
  } else {
    m_growth[ a_i ]->m_lownut = false;
  }
}
 
 
PlantGrowthData::PlantGrowthData( const char* a_vegcurvefile )
{
  // Just finds out how many veg curves there are.
  ifstream ist(a_vegcurvefile, ios::in);
  unsigned int NoPlants = FindCropNum(ist);
 
  for ( unsigned int i=0; i<NoPlants; i++) {
    CropGrowth* temp;
    temp = new CropGrowth;
    m_growth[ i ] = temp;
    SetVegNum( i, ist, a_vegcurvefile );
 
/*
    for (int mmm=0; mmm < 10; ++mmm) {
        int myin;
        ist >> myin;
        cout <<
    }
*/
 
    for (unsigned int j=0; j<5; j++) { // for each growth phase
      // 'Local' index into crop growth curves.
      int lk = 0;
      for (unsigned int k=0; k<MaxNoInflections; k++) {
        // for each inflection point
        int entry;
        //fscanf( m_ifile, "%d", &entry );
        ist >> entry;
        float f1=0,f2=0,f0=0;
        if ( entry == -1 ) {
          // Crop start data.
          m_growth[ i ]->m_start_valid[j] = true;
          //fscanf( m_ifile, "%g %g %g",&f1,&f0,&f2);
          ist >> f1 >> f0 >> f2;
                  FloatToDouble(m_growth[ i ]->m_start[j][1],f1);
                  FloatToDouble(m_growth[ i ]->m_start[j][0],f0);
                  FloatToDouble(m_growth[ i ]->m_start[j][2],f2);
        } else {
          // Add inflection point to normal growth curves.
          m_growth[ i ]->m_dds[j][lk] = (double)entry;
      //fscanf( m_ifile, "%g %g %g",&f1,&f0,&f2);
      ist >> f1 >> f0 >> f2;
          FloatToDouble(m_growth[ i ]->m_slopes[j][1][lk],f1);
          FloatToDouble(m_growth[ i ]->m_slopes[j][0][lk],f0);
          FloatToDouble(m_growth[ i ]->m_slopes[j][2][lk],f2);
          lk++;
        }
      } // MaxNoInflections
    } // Growth Phases
  } // NoPlants
  //fclose( m_ifile );
}
 
PlantGrowthData::~PlantGrowthData()
{
  for ( unsigned int i=0; i<m_growth.size(); i++ )
    delete m_growth[i];
}
 
CropGrowth::CropGrowth( void )
{
  for (unsigned int j=0; j<5; j++) {
    m_start_valid[j] = false;
    for ( unsigned int k=0; k<3; k++) {
      m_start[j][k] = 0.0;
    }
  }
}
 
int PlantGrowthData::VegTypeToCurveNum(TTypesOfVegetation VegReference)
{
  char error_num[20];
 
  switch (VegReference)
    {
    case tov_OSpringBarleyPigs:
    case tov_OSpringBarley:
    case tov_OSpringBarleyExt:
                return 101;
    case tov_SpringBarleySilage:
        case tov_SpringBarley:
        case tov_SpringBarleySpr:
        case tov_SpringBarleyPTreatment:
        case tov_SpringBarleySKManagement:
    case tov_SpringBarleyStrigling:
    case tov_SpringBarleyStriglingSingle:
    case tov_SpringBarleyStriglingCulm:
        case tov_NorwegianSpringBarley:
        case tov_PLSpringBarley:        //Needs to be changed later
        case tov_PLSpringWheat:
        case tov_PLSpringBarleySpr:
        case tov_NLSpringBarley:
        case tov_NLSpringBarleySpring:
      return 1;
    case tov_WinterBarley:
    case tov_WinterBarleyStrigling:
        case tov_PLWinterBarley:        //Needs to be changed later
      return 2;
    case tov_OWinterBarley:
    case tov_OWinterBarleyExt:
      return 102;
    case tov_WinterWheat:
    case tov_WinterWheatShort:
    case tov_WinterWheatStrigling:
    case tov_WinterWheatStriglingCulm:
    case tov_WinterWheatStriglingSingle:
    case tov_AgroChemIndustryCereal:
    case tov_WWheatPControl:
    case tov_WWheatPToxicControl:
    case tov_WWheatPTreatment:
        case tov_PLWinterWheat: // Needs to be changed later
        case tov_PLWinterWheatLate:
        case tov_DummyCropPestTesting:  // just for testing of spraying distribution
        case tov_NLWinterWheat:
      return 4;
        case tov_OWinterWheat:
        case tov_OWinterWheatUndersown:
                return 104;
    case tov_WinterRye:
    case tov_WinterRyeStrigling:
        case tov_PLWinterRye:   //Needs to be changed later
      return 5;
    case tov_OWinterRye:
    return 105;
    case tov_Oats:
        case tov_NorwegianOats:
      return 6;
    case tov_OOats:
      return 106;
    case tov_Maize:
    case tov_MaizeSilage:
    case tov_MaizeStrigling:
        case tov_PLMaize:       //Needs to be changed later
        case tov_PLMaizeSilage: //Needs to be changed later
        case tov_NLMaize:
        case tov_NLMaizeSpring:
      return 8;
    case tov_OMaizeSilage:
          return 108;
    case tov_SpringBarleyCloverGrass:
    case tov_SpringBarleyCloverGrassStrigling:
    case tov_SpringBarleySeed:
    case tov_SpringBarleyPeaCloverGrassStrigling:
      return 13;
    case tov_OBarleyPeaCloverGrass:
    case tov_OSBarleySilage:
    case tov_OSpringBarleyGrass:
    case tov_OSpringBarleyClover:
      return 113;
    case tov_WinterRape:
    case tov_WinterRapeStrigling:
        case tov_PLWinterRape: // Needs to be changed later
      return 22;
    case tov_OWinterRape:
      return 22;
        case tov_PermanentGrassLowYield:
        case tov_PermanentGrassTussocky:
                return 25;
        case tov_PermanentGrassGrazed:
        case tov_OPermanentGrassGrazed:
      return 26;
    case tov_SeedGrass1:
    case tov_SeedGrass2:
        case tov_OSeedGrass1:
        case tov_OSeedGrass2:
                return 27;
    case tov_FodderGrass:
    case tov_CloverGrassGrazed1:
    case tov_CloverGrassGrazed2:
    case tov_OCloverGrassGrazed1:
    case tov_OCloverGrassGrazed2:
    case tov_OCloverGrassSilage1:
    case tov_OGrazingPigs:
        case tov_PLFodderLucerne2:      //Needs to be changed later
        case tov_NLGrassGrazed2:
        case tov_NLGrassGrazedLast:
        case tov_NLPermanentGrassGrazed:
        return 29;
        case tov_PLFodderLucerne1:      //Needs to be changed later
        case tov_NLGrassGrazed1:
        case tov_NLGrassGrazed1Spring:
                return 70;
    case tov_OFieldPeas:
    case tov_OFieldPeasSilage:
    case tov_FieldPeas:
    case tov_FieldPeasStrigling:
        case tov_FieldPeasSilage:
        case tov_BroadBeans:
        case tov_PLBeans:
        case tov_NLCatchPeaCrop:
      return 30;
    case tov_Carrots:
        case tov_PLCarrots:     //Needs to be changed later
        case tov_NLCarrots:
        case tov_NLCabbage:     //Needs to be changed later
        case tov_NLCarrotsSpring:
        case tov_NLCabbageSpring:       //Needs to be changed later
      return 41;
    case tov_OCarrots:
      return 141;
    case tov_Potatoes:
        case tov_NorwegianPotatoes:
        case tov_PLPotatoes:    //Needs to be changed later
        case tov_PotatoesIndustry:
        case tov_NLPotatoes:
        case tov_NLPotatoesSpring:
                return 50;
        case tov_OPotatoes:
                return 150;
        case tov_SugarBeet:
        case tov_FodderBeet:
        case tov_OFodderBeet:
        case tov_PLBeet:        //Needs to be changed later
        case tov_PLBeetSpr:
        case tov_NLBeet:
        case tov_NLBeetSpring:
                return 60;
      // Special growth mode for green but unused elements.
      // tov_PermanentSetaside Does not change growth phase no matter
      // how hard one tries to do just that.
    case tov_PermanentSetaside: return 92;
        case tov_Heath:
    case tov_Setaside:
    case tov_OSetaside: return 112;
        case tov_OrchardCrop:
        case tov_YoungForest:
        case tov_NaturalGrass:
        case tov_Wasteland:
        case tov_WaterBufferZone:
                return 90;
        case tov_NoGrowth: 
        case tov_PlantNursery:
                return 91;
        case tov_Lawn: return 94;
        case tov_OTriticale:
        case tov_Triticale:
        case tov_PLWinterTriticale:     //Needs to be changed later
                return 7;
    case tov_SpringRape: return 21;
        case tov_NLTulips:      return 80;
 
 
        default: // No matching code so we need an error message of some kind
      sprintf( error_num, "%d", VegReference );
      g_msg->Warn( WARN_FILE,
                   "PlantGrowthData::VegTypeToCurveNum(): Unknown vegetation type:",
                   error_num );
      exit( 1 );
    }
}
 
 
bool  PlantGrowthData::StartValid( int a_veg_type, int a_phase )
{
  int a=m_numbers[ a_veg_type ];
  CropGrowth* p=m_growth[a];
  return p-> m_start_valid[ a_phase ];
}
 
PollenNectarDevelopmentData::PollenNectarDevelopmentData(string a_toleinputfile, string a_tovinputfile, Landscape* a_land)
{
        // This forms an empty entry for no pollen or nectar tov types
        vector<double> empty;
        vector<int> index;
        for (int i = 0; i < 366; i++) {
                empty.push_back(0.0);
                index.push_back(i);
        }
        // make space in the lookup table for all tov types
        m_tov_pollencurvetable.resize(tov_Undefined);
        for (int i = 0; i < int(tov_Undefined); i++) m_tov_pollencurvetable[i] = 0;
        // make space in the lookup table for all tole types
        m_tole_pollencurvetable.resize(tole_Foobar);
        for (int i = 0; i < int(tole_Foobar); i++) m_tole_pollencurvetable[i] = 0;
 
        if (cfg_pollen_nectar_on.value())
        {
                ifstream infile(a_tovinputfile.c_str(), ios::in);
                //check if there is an input file
                if (!infile.is_open()) {
                        g_msg->Warn("PollenNectarDevelopmentData::PollenNectarDevelopmentData Cannot open the file", a_toleinputfile.c_str());
                        exit(1);
                }
                int no_curves;
                int tov, curvenum;
                infile >> no_curves;
                for (int i = 0; i < no_curves-1; i++)
                {
                        infile >> tov >> curvenum;
                        // convert from the reference number to tov type
                        tov = a_land->TranslateVegTypes(tov);
                        // store the reference in the pollen nectar lookup table
                        m_tov_pollencurvetable[tov] = curvenum;
                }
                // Now read the curve themselves
                m_tovPollenCurves.resize(no_curves);
                m_tovNectarCurves.resize(no_curves);
                // each curve is composed to two sets of 365 numbers. The first being the index set - always 1-365 so there is no need to read this in, so it is filled in here
                vector<double> slopes(366, 0);
                PollenNectarDevelopmentCurve * p_curve;
                p_curve = new PollenNectarDevelopmentCurve(&index, &empty);
                m_tovPollenCurves[0] = p_curve;
                m_tovNectarCurves[0] = p_curve;
                for (int i = 0; i < no_curves; i++)
                {
                        // Read the curve number
                        infile >> curvenum; 
                        infile >> slopes[365];
                        // Read in the pollen data and save
                        for (int d = 0; d < 365; d++) infile >> slopes[d];
                        p_curve = new PollenNectarDevelopmentCurve(&index, &slopes);
                        m_tovPollenCurves[curvenum] = p_curve;
                        // Read in the nectar data and save
                        infile >> slopes[365];
                        for (int d = 0; d < 365; d++) infile >> slopes[d];
                        p_curve = new PollenNectarDevelopmentCurve(&index, &slopes);
                        m_tovNectarCurves[curvenum] = p_curve;
                }
                infile.close();
                infile.open(a_tovinputfile.c_str(), ios::in);
                //check if there is an input file
                if (!infile.is_open()) {
                        g_msg->Warn("PollenNectarDevelopmentData::PollenNectarDevelopmentData Cannot open the file", a_toleinputfile.c_str());
                        exit(1);
                }
                infile >> no_curves;
                int tole, tolecurvenum;
                for (int i = 0; i < no_curves-1; i++)
                {
                        infile >> tole >> tolecurvenum;
                        // convert from the reference number to tov type
                        tole = a_land->TranslateVegTypes(tole);
                        // store the reference in the pollen nectar lookup table
                        m_tole_pollencurvetable[tole] = tolecurvenum;
                }
                // Now read the curve themselves
                m_tolePollenCurves.resize(no_curves);
                m_toleNectarCurves.resize(no_curves);
                p_curve = new PollenNectarDevelopmentCurve(&index, &empty);
                m_tolePollenCurves[0] = p_curve;
                m_toleNectarCurves[0] = p_curve;
                for (int i = 1; i <= no_curves; i++)
                {
                        // Read the curve number
                        infile >> tolecurvenum;
                        // Read in the pollen data and save
                        infile >> slopes[366];
                        for (int d = 0; d < 365; d++) infile >> slopes[d];
                        p_curve = new PollenNectarDevelopmentCurve(&index, &slopes);
                        m_tolePollenCurves[tolecurvenum] = p_curve;
                        // Read in the nectar data and save
                        infile >> slopes[365];
                        for (int d = 0; d < 365; d++) infile >> slopes[d];
                        p_curve = new PollenNectarDevelopmentCurve(&index, &slopes);
                        m_toleNectarCurves[tolecurvenum] = p_curve;
                }
                infile.close();
        }
        else
        {
                PollenNectarDevelopmentCurve*  pnc = new PollenNectarDevelopmentCurve(&index, &empty);
                // Set all to zero curve
                m_tovPollenCurves.push_back(pnc);
                m_tovNectarCurves.push_back(pnc);
                m_tolePollenCurves.push_back(pnc);
                m_toleNectarCurves.push_back(pnc);
                for (int tov = 0; tov<int(tov_Undefined); tov++)
                {
                        m_tov_pollencurvetable[tov] = 0;
                }
                for (int tole = 0; tole<int(tole_Foobar); tole++)
                {
                        m_tole_pollencurvetable[tole] = 0;
                }
        }
}
 
PollenNectarDevelopmentData::~PollenNectarDevelopmentData()
{
        for (int i = 0; i < m_tovPollenCurves.max_size(); i++) delete m_tovPollenCurves[i];
        for (int i = 0; i < m_tovNectarCurves.max_size(); i++) delete m_tovNectarCurves[i];
        for (int i = 0; i < m_tolePollenCurves.max_size(); i++) delete m_tolePollenCurves[i];
        for (int i = 0; i < m_toleNectarCurves.max_size(); i++) delete m_toleNectarCurves[i];
}
 
PollenNectarDevelopmentCurveSet PollenNectarDevelopmentData::GetPollenNectarCurvePtr(int a_almassLEref)
{
        PollenNectarDevelopmentCurveSet set;
        set.m_pollencurveptr = m_tolePollenCurves[0];;
        set.m_nectarcurveptr = m_toleNectarCurves[0];;
        if ((a_almassLEref == 0) || (g_letype->TranslateEleTypes(a_almassLEref) == tole_Field) || (g_letype->TranslateEleTypes(a_almassLEref) == tole_UnsprayedFieldMargin)) return set;
        set.m_pollencurveptr = m_tolePollenCurves[m_tole_pollencurvetable[g_letype->TranslateEleTypes(a_almassLEref)]];
        set.m_nectarcurveptr = m_toleNectarCurves[m_tole_pollencurvetable[g_letype->TranslateEleTypes(a_almassLEref)]];
        return set;
}
 
PollenNectarDevelopmentCurveSet PollenNectarDevelopmentData::tovGetPollenNectarCurvePtr(int a_tov_ref)
{
        PollenNectarDevelopmentCurveSet set;
        set.m_pollencurveptr = m_tolePollenCurves[0];;
        set.m_nectarcurveptr = m_toleNectarCurves[0];;
        set.m_pollencurveptr = m_tovPollenCurves[m_tov_pollencurvetable[ a_tov_ref ]];
        set.m_nectarcurveptr = m_tovNectarCurves[m_tov_pollencurvetable[a_tov_ref]];
        return set;
}
 
PollenNectarQuality::PollenNectarQuality()
{
        m_quantity = 0;
        m_quality = 0;
 
}
 
PollenNectarQuality::PollenNectarQuality(double a_quantity, double a_quality)
{
        m_quantity = a_quantity; 
        m_quality = a_quality;
}