pesticide.cpp

Go to the documentation of this file.

//
// pesticide.cpp
//
/*
*******************************************************************************************************
Copyright (c) 2003, Christopher John Topping, EcoSol
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
//#define __WithinOrchardPesticideSim__
 
#include <iostream>     // std::cout
#include <algorithm>    // std::fill
#include <vector>       // std::vector
#include <cstdio> // Only needed for debugging purposes.
#include <math.h>
#include <fstream>
#include "ls.h"
 
using namespace std;
 
CfgBool l_pest_enable_pesticide_engine("PEST_ENABLE_PESTICIDE_ENGINE", CFG_CUSTOM, false);
CfgInt l_pest_NoPPPs("PEST_NO_OF_PPPS", CFG_CUSTOM, 1,1,10);
static CfgFloat l_pest_ai_half_life("PEST_AI_HALF_LIFE", CFG_CUSTOM, 10.0);
static CfgFloat l_pest_ai_half_life_Soil("PEST_AI_HALF_LIFE_SOIL", CFG_CUSTOM, 10.0);
static CfgFloat l_pest_ai_half_life_Veg("PEST_AI_HALF_LIFE_VEG", CFG_CUSTOM, 10.0);
static CfgFloat l_pest_diffusion_slope("PEST_DRIFT_SLOPE", CFG_CUSTOM, -0.6122);
static CfgInt   l_pest_diffusion_grid_count("PEST_DIFFUSION_GRID_COUNT",CFG_CUSTOM, 1 );
static CfgFloat l_pest_zero_threshold("PEST_ZERO_THRESHOLD",CFG_CUSTOM, 0.00001 );
 
CfgFloat cfg_pest_product_1_amount("PEST_PRODUCT_ONE_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_2_amount("PEST_PRODUCT_TWO_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_3_amount("PEST_PRODUCT_THREE_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_4_amount("PEST_PRODUCT_FOUR_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_5_amount("PEST_PRODUCT_FIVE_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_6_amount("PEST_PRODUCT_SIX_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_7_amount("PEST_PRODUCT_SEVEN_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_8_amount("PEST_PRODUCT_EIGHT_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_9_amount("PEST_PRODUCT_NINE_AMOUNT", CFG_CUSTOM, 1.0);
CfgFloat cfg_pest_product_10_amount("PEST_PRODUCT_TEN_AMOUNT", CFG_CUSTOM, 1.0);
 
CfgBool l_pest_use_application_rate("PEST_USE_APPLICATIONRATE", CFG_CUSTOM, false);
CfgFloat l_pest_daily_mort( "PEST_DAILY_MORTALITY", CFG_CUSTOM, 0.25 );
CfgFloat l_pest_daily_mort2( "PEST_DAILY_MORTALITY_TWO", CFG_CUSTOM, 0.25 );
CfgFloat l_pest_trigger_threshold1( "PEST_TRIGGER_THRESHOLD_ONE", CFG_CUSTOM, 1.0 );
CfgFloat l_pest_trigger_threshold2( "PEST_TRIGGER_THRESHOLD_TWO", CFG_CUSTOM, 1.0 );
 
CfgInt cfg_pest_productapplic_startdate("PEST_PRODUCTAPPLIC_STARTDATE", CFG_CUSTOM, -1);
CfgInt cfg_pest_productapplic_startdate2("PEST_PRODUCTAPPLIC_STARTDATE_TWO",CFG_CUSTOM,-1);
CfgInt cfg_pest_productapplic_startdate3("PEST_PRODUCTAPPLIC_STARTDATE_THREE",CFG_CUSTOM,-1);
CfgInt cfg_pest_productapplic_period("PEST_PRODUCTAPPLIC_PERIOD",CFG_CUSTOM,1);
 
CfgBool cfg_pest_residue_or_rate("PEST_RESIDUE_OR_RATE",CFG_CUSTOM,true);
 
CfgBool cfg_pest_springbarley_on("PEST_SPRINGBARLEY_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_springwheat_on("PEST_SPRINGWHEAT_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_winterwheat_on("PEST_WINTERWHEAT_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_winterbarley_on("PEST_WINTERBARLEY_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_wintertriticale_on("PEST_WINTERTRITICALE_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_winterrye_on("PEST_WINTERRYE_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_winterrape_on("PEST_WINTERRAPE_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_potatoes_on("PEST_POTATOES_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_beet_on("PEST_BEET_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_cabbage_on("PEST_CABBAGE_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_carrots_on("PEST_CARROTS_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_tulips_on("PEST_TULIPS_ON", CFG_CUSTOM, false);
CfgBool cfg_pest_SBS_ERA("PEST_SBS_ERA", CFG_CUSTOM, false);
 
class Pesticide *g_pest;
extern Landscape *g_land;
 
 
 
void Pesticide::Tick(void)
{
        if (!l_pest_enable_pesticide_engine.value())
        {
                return;
        }
        // Calculate DT50 rates
        // DT50 will change by the rate of 0.22 extra per degree different from 20
        double temp = g_weather->GetTemp();
        double Dt50 = l_pest_ai_half_life.value()* exp(0.094779*(20 - temp));
 
        m_pest_daily_decay_frac = pow(10.0, log10(0.5) / Dt50);
        // Get todays rainfall
        m_rainfallcategory = (unsigned)floor(g_land->SupplyRain() * 100 + 0.5);
        if (m_rainfallcategory >= 100) m_rainfallcategory = 99;
        // Now make it an index to the 100x100 array.
        m_rainfallcategory *= 100;
        if (!l_pest_enable_pesticide_engine.value())
        {
                for (int ppp = ppp_1; ppp < m_NoPPPs; ppp++)
                {
                        DailyQueueClear((PlantProtectionProducts) ppp);
                }
                return;
        }
        for (int ppp = ppp_1; ppp < m_NoPPPs; ppp++)
        {
                MainMapDecay((PlantProtectionProducts)ppp);
                DailyQueueProcess((PlantProtectionProducts)ppp);
                DailyQueueClear((PlantProtectionProducts)ppp);
        }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::DailyQueueClear( PlantProtectionProducts a_ppp )
{
        if (! m_daily_spray_queue[a_ppp].empty())  m_daily_spray_queue[a_ppp].resize( 0 );
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::DailyQueueAdd( LE*   a_element_sprayed, double a_amount, PlantProtectionProducts a_ppp )
{
  PesticideEvent *l_event = new PesticideEvent( a_element_sprayed, a_amount, a_ppp );
    m_daily_spray_queue[a_ppp].push_back(l_event);
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::DailyQueueProcess( PlantProtectionProducts a_ppp )
{
        if (0 == m_daily_spray_queue[a_ppp].size())
    // Event queue empty, nobody sprayed anything today.
    return;
 
  // Spraying at least one field. Force the main pesticide map
  // decay method to run tomorrow.
  m_something_to_decay[a_ppp] = true;
 
  for ( unsigned int i=0; i<m_daily_spray_queue[a_ppp].size(); i++ )
  {
          m_wind = g_land->SupplyWindDirection();
          int minx=m_daily_spray_queue[a_ppp][i]->m_sprayed_elem->GetMinX();
          int maxx=m_daily_spray_queue[a_ppp][i]->m_sprayed_elem->GetMaxX();
          int miny=m_daily_spray_queue[a_ppp][i]->m_sprayed_elem->GetMinY();
          int maxy=m_daily_spray_queue[a_ppp][i]->m_sprayed_elem->GetMaxY();
          // For cover we use the crop, and save this for later
          double cover = m_daily_spray_queue[a_ppp][i]->m_sprayed_elem->GetVegCover();
          TwinMapClear(minx, miny, maxx, maxy);
          // Add the amount in m_amount to the twin map (all squares covered with the polygon get m_amount added.
          TwinMapSpray( m_daily_spray_queue[a_ppp][i]->m_sprayed_elem, m_daily_spray_queue[a_ppp][i]->m_amount, minx, miny, maxx, maxy);
          // This adds it to the main map and if necessary sorts out the allocation between veg and soil.
          TwinMapDiffusion(minx, miny, maxx, maxy, cover, a_ppp);
  }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::TwinMapClear( int a_minx, int a_miny, int a_maxx, int a_maxy )
{
        for ( int y=a_miny; y<=a_maxy; y++ )
        {
                int t = y * m_pest_map_width;
                for ( int x=a_minx; x<=a_maxx; x++ )
                {
                        m_pest_map_twin[t + x] = 0.0;
                }
        }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::TwinMapSpray( LE* a_element_sprayed, double a_amount, int a_minx, int a_miny, int a_maxx, int a_maxy )
{
        /* Replaced with more detailed fate code for EFSA June 2014
        if (l_pest_use_application_rate.value()) {
                // We are applying a field rate. The actual residue needs to be calculated here
                double biomass = a_element_sprayed->GetVegBiomass();
                double cover = a_element_sprayed->GetVegCover();
                // rate is in mg/m and needs to be converted to mg/kg
                double residue = (a_amount/(biomass*0.001)) * 0.435855; // g to kg veg  // 0.435855 is a specific calculation to obtain the same residue that Joe Crocker used for the scale of use study
                a_amount = residue * cover;
        }
        */
        double l_fractional_amount = a_amount * m_prop; //  m_prop is a constant related to the grid area (i.e. how many squares are in one grid square).
        int l_large_map_index = a_element_sprayed->GetMapIndex();
        for ( int y=a_miny; y<=a_maxy; y++ ) {
                for ( int x=a_minx; x<=a_maxx; x++ ) {
                        if ( m_land->Get( x, y ) == l_large_map_index )
                        {
                                // This adds the l_fractional_amount to the twin map
                                TwinMapSprayPixel( x, y, l_fractional_amount );
                        }
                }
        }
        TwinMapSprayCorrectBorders();
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::TwinMapSprayCorrectBorders( void )
{
  if ( m_x_excess ) {
   for ( unsigned int i=0; i<m_pest_map_width-1; i++ ) {
      m_pest_map_twin[ i * m_pest_map_width + m_pest_map_height-1 ] *=  m_corr_x;
    }
  }
 
  if ( m_y_excess ) {
    unsigned int l_additive = (m_pest_map_height-1)*m_pest_map_width;
    for ( unsigned int i=0; i<m_pest_map_height-1; i++ ) {
      m_pest_map_twin[ i + l_additive ] *=      m_corr_y;
    }
  }
 
  if ( m_x_excess && m_y_excess ) {
    m_pest_map_twin[ m_pest_map_size-1 ] *= (((double)(PEST_GRIDAREA)) / ((double)(m_x_excess*m_y_excess)));
  }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
/*
void Pesticide::TwinMapCopyToMask( int a_minx, int a_miny, int a_maxx, int a_maxy  )
{
        // Copy the newly calculated pesticide amounts over to the mask map
        for ( int y=a_miny; y<=a_maxy; y++ )
        {
                int t = y * m_pest_map_width;
                for ( int x=a_minx; x<=a_maxx; x++ )
                {
                        m_pest_map_mask[t + x]  = m_pest_map_twin[t + x];
                }
        }
}
*/
//-----------------------------------------------------------------------------------------------------------------------------
 
/*
void Pesticide::TwinMapAddToMain( void )
{
  for ( unsigned int i=0; i<m_pest_map_size; i++ ) {
    m_pest_map_main[i] += m_pest_map_twin[i];
  }
}
//-----------------------------------------------------------------------------------------------------------------------------
*/
void Pesticide::TwinMapDiffusion( int a_minx, int a_miny, int a_maxx, int a_maxy, double a_cover, PlantProtectionProducts a_ppp )
{
        for ( int y=a_miny; y<a_maxy; y++ )
        {
                int t = y*m_pest_map_width;
                for ( int x=a_minx; x<a_maxx; x++ )
                {
                        if ( m_pest_map_twin[ t + x ] > 0.0 )
                        {
                                double l_amount = m_pest_map_twin[ t + x ];
                                for ( unsigned int i=0; i<m_diffusion_mask[m_wind].size(); i++ )
                                {
                                        // This adds the pesticide to the main map - simply adds the calculated amount.
                                        DiffusionSprayPixel( x + m_diffusion_mask[m_wind][i]->Getdx(),m_pest_map_width,y + m_diffusion_mask[m_wind][i]->Getdy(),
                                                m_pest_map_height,m_diffusion_mask[m_wind][i]->GetFraction() *l_amount, a_cover, a_ppp);
                                }
                        }
                }
        }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
#ifdef __DETAILED_PESTICIDE_FATE
inline void Pesticide::DiffusionSprayPixel(int a_x, int a_limit_x, int a_y, int a_limit_y, double a_amount, double a_cover, PlantProtectionProducts a_ppp)
#else
inline void Pesticide::DiffusionSprayPixel(int a_x, int a_limit_x, int a_y, int a_limit_y, double a_amount, double /* a_cover */, PlantProtectionProducts a_ppp)
#endif
{
        // First we have to do the typical out of bounds checks - if these fail do nothing, the pesticide fell off the world
        if (a_x < 0 || a_x >= a_limit_x || a_y < 0 || a_y >= a_limit_y) return;
        // Now calculate the coordinate entry in the array and store this as l_coord
        int l_coord = a_y * a_limit_x + a_x;
#ifdef __DETAILED_PESTICIDE_FATE
        // Here we need to calculate the partition of pesticide into two compartments.
        m_pest_map_vegcanopy[a_ppp][l_coord] += a_amount * a_cover;
        m_pest_map_soil[a_ppp][l_coord] += a_amount - m_pest_map_vegcanopy[a_ppp][l_coord];
#else
        m_pest_map_main[a_ppp][l_coord] += a_amount;
#endif
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::DiffusionMaskInit( void )
{
  int l_grid        = l_pest_diffusion_grid_count.value();
  const int l_side_length = l_grid * 2 + 1;
  double* l_diff;
  l_diff = new double[l_side_length];
  int cells = 1;
  double sum = 1.0;
  l_diff[0] = 1.0;
  for (int l=1; l<l_side_length; l++)
  {
        l_diff[l] = DiffusionFunction(l);
        sum += l_diff[l];
  }
  for (int l=0; l<l_side_length; l++)
  {
        l_diff[l] /= sum;
        l_diff[l] /= cells;
        cells += 2; // 2 more cells per row away from centre square
  }
  // Calculated for 4 wind directions
  for (int wind=0; wind < 4; wind ++)
  {
          m_diffusion_mask[wind].resize( l_side_length * l_side_length );
          for ( int x = 0; x< l_side_length; x++ )
          {
                for ( int y= 0; y< l_side_length; y++ )
                {
                  m_diffusion_mask[wind][x + y * l_side_length] = new Diffusor( x, y, 0);  // first zero all values
                }
          }
  }
  int strtx = l_grid;
  int strty = l_grid;
  int fin = 1;
  // North
  for (int step = 0; step <= l_grid; step++)
     {
          for (int cc=0; cc<fin; cc++)
          {
                  m_diffusion_mask[0][ strtx + (cc-step) + ((strty-step)*l_side_length)]->SetFraction(l_diff[step]);
          }
  fin += 2;
  }
  // South
  fin = 1;
  for (int step = 0; step <= l_grid; step++)
     {
          for (int cc=0; cc<fin; cc++)
          {
                  m_diffusion_mask[2][ strtx + (cc-step) + ((strty+step)*l_side_length)]->SetFraction(l_diff[step]);
          }
          fin += 2;
  }
  // East
  fin = 1;
  for (int step = 0; step <= l_grid; step++)
     {
          for (int cc=0; cc<fin; cc++)
          {
                  m_diffusion_mask[3][ (strtx + step) + ((strty+(cc-step))*l_side_length)]->SetFraction(l_diff[step]);
          }
          fin += 2;
  }
  // West
  fin = 1;
  for (int step = 0; step <= l_grid; step++)
     {
          for (int cc=0; cc<fin; cc++)
          {
                  m_diffusion_mask[1][ (strtx - step) + ((strty+(cc-step))*l_side_length)]->SetFraction(l_diff[step]);
          }
          fin += 2;
  }
 
  delete[] l_diff;
}
//-----------------------------------------------------------------------------------------------------------------------------
 
 
 
/*
void Pesticide::DiffusionMaskInit( void )
{
  int l_grid        = l_pest_diffusion_grid_count.value();
  int l_side_length = l_grid * 2 + 1;
 
  m_diffusion_mask.resize( l_side_length * l_side_length );
 
  unsigned int i    = 0;
  //double l_diff_zero = DiffusionFunction( 0.0 );
  double l_diff_sum = 0.0;
  for ( int x = -l_grid; x < l_grid + 1; x++ ) {
    for ( int y = -l_grid; y < l_grid + 1; y++ ) {
      double l_dist = sqrt( (double) x*x + y*y ); //  ***CJT*** added 28/8/2009
          // m_diffusion_mask has the square to drift spray to with the proportions relative to the amount received in the centre square
          // but the total may be greater than one, so we need the adjustment provided by the second loop
          double l_diff = DiffusionFunction(l_dist);
          l_diff_sum += l_diff;
      m_diffusion_mask[i] = new Diffusor( x, y, l_diff);  // ***CJT*** changed 28/8/2009 l_diff_zero
      i++;
    }
  }
  for ( unsigned x = 0; x < i; x++) {
                m_diffusion_mask[x]->m_fraction /= l_diff_sum;
  }
}
 
*/
 
double Pesticide::DiffusionFunction( double a_dist_meters )
{
  double pp;
//   y=(2.753767)*(x+(1.86976))**(-2.121563)
#ifdef __WithinOrchardPesticideSim__
  if (a_dist_meters==0) pp=0.7784;
    else pp=0.0277;
#else
        //pp = (2.7705 * pow(a_dist_meters, -0.9787)) * 0.01; // german drif equation (reference)
        pp = (2.7705 * pow(a_dist_meters, -1.7)) * 0.01; // german drif equation - reduction of drift by 50%
        //pp = (2.7705 * pow(a_dist_meters, -2.5)) * 0.01; // german drif equation - reduction of drift by 75%
        //pp = (2.7705 * pow(a_dist_meters, -3.5)) * 0.01; // german drif equation - reduction of drift by 90%
 
    //pp=2.753767 * pow(( a_dist_meters+1.86976 ), -2.121563); // based on drift data from the spray drift calculator within FOCUS's surface water scenarios SWASH software (FOCUS, 2001)
    //pp= exp( l_pest_diffusion_slope.value() * a_dist_meters );
    pp=pow(( a_dist_meters+1 ), l_pest_diffusion_slope.value() );
    //pp = (2.7593 * pow(a_dist_meters,-0.9778)) * 0.001;
    if (pp<(l_pest_zero_threshold.value()/10)) pp=0;  // Don't bother with the little ones
#endif
   return pp;
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::MainMapDecay(PlantProtectionProducts a_ppp)
{
        if ( ! m_something_to_decay[a_ppp] )    return;
        double l_zero         = l_pest_zero_threshold.value();
        m_something_to_decay[a_ppp] = false;
        for ( unsigned int i=0; i<m_pest_map_size; i++ )
        {
#ifdef __DETAILED_PESTICIDE_FATE
                if ((m_pest_map_vegcanopy[a_ppp][i] > l_zero) )
                {
                        //Calculate wash-off  m_pest_map_width * m_pest_map_height y * m_pest_map_width + x
                        int x = ( i % m_pest_map_width ) << PEST_GRIDSIZE_POW2;
                        int y = ( i / m_pest_map_width ) << PEST_GRIDSIZE_POW2;
                        double cover = g_landscape_p->SupplyVegCover(x,y);
                        unsigned cov = 100 * cover; // Cover is steps of zero to 99, may need to check for bounds here (may get cover of 1.0 ?)
                        double Rwp = m_pest_map_vegcanopy[a_ppp][i] *= m_RainWashoffFactor[m_rainfallcategory+cov]; // m_RainWashoffFactor is the index to the array
                        m_pest_map_soil[a_ppp][i] = (m_pest_map_soil[a_ppp][i] + Rwp) * m_pest_daily_decay_frac_Soil;
                        m_pest_map_vegcanopy[a_ppp][i] -= Rwp;
                        m_pest_map_vegcanopy[a_ppp][i] *= m_pest_daily_decay_frac_Veg;
                        m_something_to_decay[a_ppp] = true;
                }
                else
                {
                        if (m_pest_map_soil[a_ppp][i] > l_zero)
                        {
                                m_pest_map_soil[a_ppp][i] = m_pest_map_soil[a_ppp][i] * m_pest_daily_decay_frac_Soil;
                                m_something_to_decay[a_ppp] = true;
                        }
                        else
                        {
                                m_pest_map_vegcanopy[a_ppp][i] = 0.0;
                                m_pest_map_soil[a_ppp][i] = 0.0;
                        }
                }
#else
          if (m_pest_map_main[a_ppp][i] > l_zero)
          {
                  m_pest_map_main[a_ppp][i]       *= m_pest_daily_decay_frac;
                  m_something_to_decay[a_ppp]      = true;
          }
          else
          {
                  m_pest_map_main[a_ppp][i]        = 0.0;
          }
#endif
  }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
Pesticide::Pesticide( RasterMap* a_land, Landscape *a_map )
{
        if ( l_pest_enable_pesticide_engine.value())
        {
                m_NoPPPs = l_pest_NoPPPs.value();
                m_map = a_map;
                m_land = a_land;
                for (int i=0; i< m_NoPPPs; i++) m_something_to_decay[i] = false;
                m_daily_spray_queue.resize(m_NoPPPs);
 
                // Figure out critical border coordinates,
                // proportional fractions etc.
                m_prop     = 1.0 / ((double)(PEST_GRIDAREA));
                m_x_excess = m_land->MapWidth() & (PEST_GRIDSIZE-1);
                if ( m_x_excess ) {
                        m_corr_x   = (double)PEST_GRIDSIZE /(double)m_x_excess;
                }
                m_y_excess = m_land->MapHeight() & (PEST_GRIDSIZE-1);
                if ( m_y_excess ) {
                        m_corr_y = (double)PEST_GRIDSIZE / (double)m_y_excess;
                }
 
                m_pest_map_width = m_land->MapWidth() >> PEST_GRIDSIZE_POW2; // landscape width divided by the sqrt(cellsize)  e.g. cellsize 4 and landscape width = 10000, result is 10000/(sqrt(4) = 5000
                if ( m_land->MapWidth() & (PEST_GRIDSIZE-1))
                        m_pest_map_width++;
 
                m_pest_map_height = m_land->MapHeight() >> PEST_GRIDSIZE_POW2;
                if ( m_land->MapHeight() & (PEST_GRIDSIZE-1))
                        m_pest_map_height++;
 
                m_pest_map_size = m_pest_map_width * m_pest_map_height;
 
                m_pest_map_twin = new double[m_pest_map_size];
#ifdef __DETAILED_PESTICIDE_FATE
                m_pest_map_vegcanopy.resize(m_NoPPPs);
                m_pest_map_soil.resize(m_NoPPPs);
#else
                m_pest_map_main.resize(m_NoPPPs);
#endif
                for (int p = 0; p < m_NoPPPs; p++) {
#ifdef __DETAILED_PESTICIDE_FATE
                        m_pest_map_vegcanopy[p].resize(m_pest_map_size);
                        m_pest_map_soil[p].resize(m_pest_map_size);
                        fill(m_pest_map_vegcanopy[p].begin(), m_pest_map_vegcanopy[p].end(), 0.0);
                        fill(m_pest_map_soil[p].begin(), m_pest_map_soil[p].end(), 0.0);
#else
                        m_pest_map_main[p].resize(m_pest_map_size);
                        fill(m_pest_map_main[p].begin(), m_pest_map_main[p].end(), 0.0);
                
#endif
                        m_daily_spray_queue[p].resize(0);
                }
                
                for (unsigned int i = 0; i < m_pest_map_size; i++) m_pest_map_twin[i] = 0.0;    
 
                m_pest_daily_decay_frac = pow(10.0, log10(0.5) / l_pest_ai_half_life.value());
#ifdef __BORLANDC__
                m_pest_daily_decay_frac_Veg = log(double(2.0))/l_pest_ai_half_life_Veg.value();
#else
                m_pest_daily_decay_frac_Veg = log(2.0)/l_pest_ai_half_life_Veg.value();
#endif
 
                m_pest_daily_decay_frac_Soil = pow(10.0, log10(0.5) / l_pest_ai_half_life_Soil.value());
                // This calculates
                CalcRainWashOffFactors();
                DiffusionMaskInit();
                DiffusionMaskInitTest();
        }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
void Pesticide::CalcRainWashOffFactors()
{
        double rainsteps = 0.1; // mm, multiply rainfall by 100 to get this, save result as integer
        double coversteps = 1.0; // %, multiply cover by 100 to get this, save result as integer
        for (int i = 0; i < 100; i++)
        {
                double SC = i * coversteps;
                double LAI = log(0 - (1 - SC))*1.666666667; // This is the inverse of Beer's Law
                for (int r = 0; r < 100; r++)
                {
                        double P = r * rainsteps;
                        double Pi = LAI*(1 - (1 / ((1 + SC * P) / LAI)));
                        double Rw = 0.25*(SC * P - Pi); // Rw is a proportion of canopy pesticide washed off.
                        m_RainWashoffFactor[r * 100 + i] = Rw;
                }
        }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
Pesticide::~Pesticide(void)
{
        if (l_pest_enable_pesticide_engine.value())
        {
                for (int winddir = 0; winddir < 4; winddir++)
                {
                        for (unsigned int i = 0; i < m_diffusion_mask[winddir].size(); i++)
                        {
                                delete m_diffusion_mask[winddir][i];
                        }
                        m_diffusion_mask[winddir].resize(0);
                }
                delete m_pest_map_twin;
        }
}
//-----------------------------------------------------------------------------------------------------------------------------
 
#ifdef PEST_DEBUG
 
void Pesticide::DiffusionMaskInitTest( void )
{
  int l_grid        = l_pest_diffusion_grid_count.value();
  const int l_side_length = l_grid * 2 + 1;
 
 // Need a debug test, output all diffusion masks to file
  ofstream ofile("diffusionmask.txt",ios::out);
  for (int i=0; i<4; i++)
  {
          for ( int x = 0; x< l_side_length; x++ )
          {
                  for ( int y= 0; y< l_side_length; y++ )
                  {
                          ofile << m_diffusion_mask[i][x+(y*l_side_length)]->GetFraction() << '\t';
                  }
                  ofile << endl;
          }
          ofile << endl;
          ofile << endl;
  }
  ofile.close();
 }
//-----------------------------------------------------------------------------------------------------------------------------
 
// For file saving and loading.
#define CFG_CHANNEL_BITS      8
#define CFG_CHANNEL_MAXVAL    (2^CFG_CHANNEL_BITS-1)
 
#define SV_UINT32 unsigned int
#define SV_INT32  int
#define SV_UINT8  unsigned char
#define SV_INT8   char
 
bool Pesticide::SavePPM( double *a_map,
                         int a_beginx, int a_width,
                         int a_beginy, int a_height,
                         char* a_filename )
{
  a_beginx = 0;
  a_width  = m_pest_map_width;
  a_beginy = 0;
  a_height = m_pest_map_height;
 
  SV_UINT32 linesize   = a_width*3;
  SV_UINT8* linebuffer = (SV_UINT8*)malloc(sizeof(SV_UINT8)* linesize);
 
  if ( linebuffer == NULL ) {
    g_msg->Warn( WARN_FATAL,
                 "Pesticide::SavePPM(): Out of memory!", "" );
    exit(1);
  }
 
  FILE* l_file;
  l_file=fopen(a_filename, "w" );
  if ( !l_file ) {
    printf("PesticideTest::SavePPM(): "
           "Unable to open file for writing: %s\n",
           a_filename );
    exit(1);
  }
 
  fprintf( l_file, "P6\n%d %d %d\n",
           a_width,
           a_height,
           255 );
 
  for ( int line=a_beginy; line< a_beginy + a_height; line++ ) {
    int i = 0;
    for ( int column=a_beginx; column < a_beginx + a_width; column++ ) {
      int localcolor = (int)( a_map[ line * m_pest_map_width + column ]
                                  * 255.0);
      if ( localcolor <= 255 ) {
        linebuffer [ i++ ] =  char (localcolor        & 0xff);
        linebuffer [ i++ ] = 0;
        linebuffer [ i++ ] = 0;
      } else {
        linebuffer [ i++ ] = 255;
        localcolor -= 255;
        if ( localcolor <= 255 ) {
          linebuffer [ i++ ] = char (localcolor);
          linebuffer [ i++ ] = 0;
        } else {
          linebuffer [ i++ ] = 255;
          localcolor -= 255;
          if ( localcolor <= 255 ) {
            linebuffer [ i++ ] = char (localcolor);
          } else {
            linebuffer [ i++ ] = 255;
          }
        }
      }
    }
    fwrite( linebuffer, sizeof(SV_UINT8), linesize, l_file );
  }
 
  fclose( l_file );
  free( linebuffer );
  return true;
}
//-----------------------------------------------------------------------------------------------------------------------------
 
/*
void Pesticide::Test( Landscape *a_map )
{
  for ( int i=0; i<730; i++) {
 
    if ( random(100) < 10 ) {
      int l_sprays = random(5)+1;
      for ( int j=0; j<l_sprays; j++ ) {
        int l_x = random(1500);
        int l_y = random(1500);
        DailyQueueAdd( a_map->
                       SupplyLEPointer( a_map->SupplyPolyRef( l_x, l_y )),
                       1.0 );
      }
      char l_filename[20];
      sprintf( l_filename, "ppms/p%04d.ppm", i );
      SavePPM( m_pest_map_main,
               0, 0, 0, 0,
               l_filename );
    }
 
    Tick();
  }
}
//-----------------------------------------------------------------------------------------------------------------------------
*/
#endif // PEST_DEBUG
 
PesticideMap::PesticideMap(int a_startyear, int a_noyears, int a_cellsize, Landscape * a_landscape, RasterMap* a_land, bool a_typeofmap)
{
        m_typeofmap = a_typeofmap;
        m_startyear = a_startyear;
        m_endyear = a_startyear + a_noyears;
        m_cellsize = a_cellsize;
        m_OurLandscape = a_landscape;
        m_Rastermap = a_land;
        m_pmap_width = m_OurLandscape->SupplySimAreaWidth() / m_cellsize;
        m_pmap_height = m_OurLandscape->SupplySimAreaHeight() / m_cellsize; ;
 
        // Need to create the grid to fill. This grid is always used.
        m_pmap_insecticides = new vector<double>;
        m_pmap_insecticides->resize(m_pmap_width*m_pmap_height);
        fill(m_pmap_insecticides->begin(), m_pmap_insecticides->end(), 0);
        ofstream * m_PMap;
        m_PMap = new ofstream("ALMaSS_InsecticideMap.txt", ios::out);
        (*m_PMap) << "Years" << ' ' << a_noyears << ' ' << "Height" << ' ' << m_pmap_height << ' ' << "Width" << ' ' << m_pmap_width << ' ' << endl;
        m_PMap->close();
        // These two grids are used only if we want all three types recoreded.
        if (!m_typeofmap)
        {
                // We need all three in this case
                m_pmap_fungicides = new vector<double>;
                m_pmap_fungicides->resize(m_pmap_width*m_pmap_height);
                fill(m_pmap_fungicides->begin(), m_pmap_fungicides->end(), 0);
                m_PMap = new ofstream("ALMaSS_FungicideMap.txt", ios::out);
                (*m_PMap) << "Years" << ' ' << a_noyears << ' ' << "Height" << ' ' << m_pmap_height << ' ' << "Width" << ' ' << m_pmap_width << ' ' << endl;
                m_PMap->close();
                m_pmap_herbicides = new vector<double>;
                m_pmap_herbicides->resize(m_pmap_width*m_pmap_height);
                fill(m_pmap_herbicides->begin(), m_pmap_herbicides->end(), 0);
                m_PMap = new ofstream("ALMaSS_HerbicideMap.txt", ios::out);
                (*m_PMap) << "Years" << ' ' << a_noyears << ' ' << "Height" << ' ' << m_pmap_height << ' ' << "Width" << ' ' << m_pmap_width << ' ' << endl;
                m_PMap->close();
        }
        else
        {
                m_pmap_fungicides = NULL;
                m_pmap_herbicides = NULL;
        }
}
 
PesticideMap::~PesticideMap()
{
        delete m_pmap_insecticides;
        if (!m_typeofmap)
        {
                // We need all three in this case
                delete m_pmap_fungicides;
                delete m_pmap_herbicides;
 
        }
}
 
bool PesticideMap::DumpPMap(vector<double>* a_map)
{
        string fname;
        if (a_map == m_pmap_insecticides) fname = "ALMaSS_InsecticideMap.txt";
        else if (a_map == m_pmap_fungicides) fname = "ALMaSS_FungicideMap.txt";
        else  fname = "ALMaSS_HerbicideMap.txt";
        {
                ofstream* m_PMap = new ofstream(fname.c_str(), ios::app);
                if (!(*m_PMap).is_open())
                {
                        g_msg->Warn(WARN_FILE, "PesticideMap::DumpMap Unable to open file for append: ", fname.c_str());
                        exit(1);
                }
                // File is OK, save the map
                (*m_PMap) << "Year " << g_date->GetYear() << "Month " << g_date->GetMonth() << endl;
                for (int y = 0; y < m_pmap_height; y++)
                {
                        for (int x = 0; x < m_pmap_width; x++)
                        {
                                (*m_PMap) << (*a_map)[y*m_pmap_width + x] << ' ';
                        }
                }
                (*m_PMap) << endl;
                m_PMap->close();
        }
        fill(a_map->begin(), a_map->end(), 0);
        return true;
}
 
 
void PesticideMap::Spray(LE* a_element_sprayed, TTypesOfPesticideCategory a_type)
{
        vector<double>* ourmap = m_pmap_insecticides; ;
        if (a_type == fungicide) ourmap = m_pmap_fungicides;
        else if (a_type == herbicide) ourmap = m_pmap_herbicides;
        int l_large_map_index = a_element_sprayed->GetMapIndex();
        int minx = a_element_sprayed->GetMinX();
        int miny = a_element_sprayed->GetMinY();
        int maxx = a_element_sprayed->GetMaxX();
        int maxy = a_element_sprayed->GetMaxY();
        for (int y = miny; y <= maxy; y++) {
                for (int x = minx; x <= maxx; x++) {
                        if (m_Rastermap->Get(x, y) == l_large_map_index)
                        {
                                // Add one m spray value of 1.0 per m
                                (*ourmap)[(x/m_cellsize)+(y/m_cellsize)*m_pmap_width]++;
                        }
                }
        }
}
//-----------------------------------------------------------------------------------------------------------------------------