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 <math.h>
#include <stdlib.h>
#include "configurator.h"
#include "maperrormsg.h"
#include "plants.h"
#include "tov_declaration.h"

extern void FloatToDouble(double&,float);
class CropData *g_crops;

double CropData::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;
}


double CropData::GetLAgreenDiff( double a_ddegs, double a_yddegs,
                int a_plant, int a_phase )
{
  return FindDiff( a_ddegs, a_yddegs, a_plant, a_phase, 0 );
}


double CropData::GetLAtotalDiff( double a_ddegs, double a_yddegs,
                int a_plant, int a_phase )
{
  return FindDiff( a_ddegs, a_yddegs, a_plant, a_phase, 1 );
}


double CropData::GetHeightDiff( double a_ddegs, double a_yddegs,
                   int a_plant, int a_phase )
{
  return FindDiff( a_ddegs , a_yddegs, a_plant, a_phase, 2 );
}



unsigned int CropData::FindCropNum( const char* a_cropcurvefile )
{
  int NoPlants;

  m_numbers.resize(201);
  for ( unsigned int i=0; i<201; i++) {
    m_numbers[ i ] = -1;
  }

  m_ifile = fopen(a_cropcurvefile, "r" );
  if (!m_ifile){
    g_msg->Warn(WARN_FILE, "CropData::CropData: Unable to open file",
        a_cropcurvefile );
    exit(1);
  }

  fscanf( m_ifile , "%d", &NoPlants );  // How many tables to read in
  m_growth.resize( NoPlants );
  m_num_crops = NoPlants;

  return NoPlants;
}



void CropData::SetVegNum( unsigned int a_i, const char* a_cropcurvefile )
{
  int ThisPlant;

  // Find out what crop and what nutrient status
  fscanf( m_ifile, "%d", &ThisPlant);

  // Check if valid plant number (from the file).
  if ( ThisPlant < 0 || ThisPlant > 200 ) {
    g_msg->Warn(WARN_FILE, "CropData::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;
  }
}


CropData::CropData( const char* a_vegcurvefile )
{
  // Just finds out how many veg curves there are.
  unsigned int NoPlants = FindCropNum( a_vegcurvefile );

  for ( unsigned int i=0; i<NoPlants; i++) {
    CropGrowth* temp;
    temp = new CropGrowth;
    m_growth[ i ] = temp;
    SetVegNum( i, a_vegcurvefile );

    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 );
    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);
          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);
      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 );
}

CropData::~CropData()
{
  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 CropData::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:
      return 1;
    case tov_WinterBarley:
    case tov_WinterBarleyStrigling:
      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:
      return 4;
    case tov_OWinterWheat:
    case tov_OWinterWheatUndersown:
        return 104;
    case tov_WinterRye:
    case tov_WinterRyeStrigling:
      return 5;
    case tov_OWinterRye:
    return 105;
    case tov_Oats:
      return 6;
    case tov_OOats:
      return 106;
    case tov_Maize:
    case tov_MaizeSilage:
    case tov_MaizeStrigling:
      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:
      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:
        return 29;
    case tov_OFieldPeas:
    case tov_OFieldPeasSilage:
    case tov_FieldPeas:
    case tov_FieldPeasStrigling:
    case tov_FieldPeasSilage:
    case tov_BroadBeans:
      return 30;
    case tov_Carrots:
      return 41;
    case tov_OCarrots:
      return 141;
    case tov_Potatoes: return 50;
    case tov_OPotatoes: return 150;
    case tov_PotatoesIndustry: return 50;
    case tov_SugarBeet:
    case tov_FodderBeet:
    case tov_OFodderBeet:
        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: return 90;
    case tov_NoGrowth: return 91;
    case tov_Lawn: return 94;
      //case tov_SpringWheat: return 3;
    case tov_OTriticale:
    case tov_Triticale: return 7;
    case tov_SpringRape: return 21;

    default: // No matching code so we need an error message of some kind
      sprintf( error_num, "%d", VegReference );
      g_msg->Warn( WARN_FILE,
           "CropData::VegTypeToCurveNum(): Unknown vegetation type:",
           error_num );
      exit( 1 );
    }
}


bool  CropData::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 ];
}