ReGen.cpp

/*
 *  ReGen.cpp
 *  ReGen 1.0.1
 *
 *  Created by Martin Koechy on Thu Feb 06 2003.
 *
 */
#include "ReGen.h"
#include <math.h>

RAINPARAMETERS::RAINPARAMETERS (int MAP, float dDMR)
{
  calcParameters (MAP, dDMR);
}

RAINPARAMETERS::RAINPARAMETERS (float A, float L, float W, float a, float l, float w, float dDMR)
{
  lA=A; lL=L; lW=W; vA=a; vL=l; vW=w; d=dDMR;
}

void RAINPARAMETERS::calcParameters (float MAP, float d) 
{ 
  // Parameter values based on monthly units
  
  // rain occurance:
  lA = 0.130034 + 0.000408 * MAP;
  lA = lA*(1-d +1.33*d*d-(0.61+1.57*lA)*d*d*d);
  lL = 5.735; // 5.235 + 0.5 months
  lW = 1.704746 + 0.000242 * MAP;
  
  // rain volume:
  vA = -13.7178 + 4.0508301 * log(MAP);
  vA = vA * (1+d);
  vL = 5.586; // 5.086 + 0.5 months
  vW = 11.336;
}


//------------------------------------------------------------------
ReGen::ReGen (RAINPARAMETERS* p_RainP)
{   
  pP = p_RainP;
  theAnnualRain = 0.0;
  
  float days[days_in_year];
  for (int d = 0; d < days_in_year; days[d] = d++);
  valarray <float> localDays (days, days_in_year);
  theDays.resize(days_in_year);
  theDays = localDays;
  
  rainProb.resize(days_in_year);

  meanRainVolume.resize(days_in_year);
  meanRainVolume = calcGauss(pP->vA, pP->vL, pP->vW, 4);

  rainyDay.resize(days_in_year);
  
  rainVolume.resize(days_in_year);
}

//------------------------------------------------------------------
void ReGen::drawYear (void)
{
  calcRain(); 
  theAnnualRain = rainVolume.sum(); // calculate right away, reuse later
}

//------------------------------------------------------------------
valarray <float> ReGen::calcGauss (float A, float L, float W, int s)
{
  // A: amplitude, L: location, W: width, s: shape
  valarray <float> G (days_in_year);
  G = theDays;
  
  G *= (12.0F/365.0F);
  G += (6.0F - L);
  for (int i = 0; i < days_in_year; i++) //  .apply does not work
        G[i] = fmod(G[i],12.0F);
  G += (L - 6.0F);
  
  G = A * exp(-pow(G - L, float(s))/float(2.0*W*W));

  return G;
}

//------------------------------------------------------------------
void ReGen::calcRain (void)
{   
  // calculate rainy days
  rainProb = calcGauss(pP->lA, pP->lL, pP->lW);

  rainyDay=0.0;
  for (int i = 0; i < days_in_year; i++)
        if (Zzg.rand_halfclosed01() < rainProb[i] && rainProb[i] > 0.05)
          rainyDay[i] = 1.0;
 
  
  // calculate rain volume on a specific day
  for (int i = 0; i < days_in_year; i++)
        rainVolume[i] = Zzg.rand_halfclosed01();
  rainVolume *= -1.0;
  rainVolume += 1.0;
  rainVolume = -log(rainVolume);
  
  rainVolume *= meanRainVolume;
  rainVolume *= rainyDay;
}

//------------------------------------------------------------------
float ReGen::getRain (int day) const
{
  return rainVolume[day];
}

//------------------------------------------------------------------
float ReGen::getAnnualRain (void) const
{
  return theAnnualRain;
}

//------------------------------------------------------------------

float ReGen::getPeriodRain (int start, int end) const
{
  float sum = 0.0;
  for (int day = start; day <= end; day++)
        sum+=rainVolume[day];
  return sum;
}

---