// Infiltration.cpp



#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <cstring>
#include <fstream>
#include <cassert>
#include <cmath>
#include "infiltration.h"

#define TEST
#undef TEST

// public functions  /////////////////////////////////////////////////////////////////////
Infiltration::Infiltration(): PowerLaw()
{
  FinalRate = 0;
  MaxStorage = 1;
  FieldCapStorage = 0;
  FinalRate = 0;
}

Infiltration::Infiltration(  float store, float coef, float exp, float max_store, 
                             float fc_store, float final  ): PowerLaw( store, coef, exp )
{
  assert( max_store >= store );
  assert( fc_store > 0 && fc_store < max_store );

  MaxStorage = max_store;
  FieldCapStorage = fc_store;
  FinalRate = final;
}


void Infiltration::get_parameters(  float store, float coef, float exp, float max_store, 
                               float fc_store, float final  )
{

  assert( store >= 0.0 );
  assert( max_store >= store );
  assert( fc_store > 0 && fc_store < max_store );

  Storage = store;
  Coeff = coef;
  MaxStorage = max_store;
  FieldCapStorage = fc_store;
  FinalRate = final;
}


float Infiltration::CalcRunoff( float InputVol, float Dt, float *DeepPerc )     
{
  float InFlux = InputVol/Dt;
  float Runoff = 0.0;
  float NegDeltaStorage = 0.0;
  float oldStorage = Storage;
  *DeepPerc = 0.0;

  if( InputVol == 0 ) {  // no runoff

    *DeepPerc = NoInput( Dt );
    return( 0.0 );

  }

  // Input > 0
  if( Storage > FieldCapStorage ) {  // must bring soil to field capacity
    Dt = FillSoilToFC( InFlux, Dt );
    if( Dt == 0.0 ) return( 0.0 );  // all input went to soil; no runoff or deep perc
  }

  // soil >= field capacity, check for ponding
  float infiltration_rate = Coeff * pow( Storage, Exponent ) + FinalRate;

  if( infiltration_rate >= InFlux ) { // no ponding yet

    if( InFlux <= FinalRate ) {  // soil losing water, will not pond
      Dt = DrainSoilToTarget( FieldCapStorage, InFlux, Dt, DeepPerc );
      if( Dt > 0 ) Storage -= InFlux  * Dt;
      return( 0.0 );
    }
    else {  // soil gaining water, may become ponded
      float ponding_volume = pow( (InFlux - FinalRate) / Coeff, 1 / Exponent );
      Dt = DrainSoilToTarget( ponding_volume, InFlux, Dt, DeepPerc );
      if( Dt == 0 ) return( 0.0 );   // no ponding or runoff
    }

   }

    // ponded conditions
    // with InputVol = 0, output returns -DeltaStorage
    NegDeltaStorage = output( 0.0, Dt );

    Runoff = (InFlux - FinalRate) * Dt - NegDeltaStorage;  // MW

    if( Runoff < 0 ) {
      Runoff = 0.0;
      Storage = (InFlux - FinalRate) * Dt + oldStorage;
    }

    *DeepPerc += FinalRate * Dt;


    //assert( Runoff >= 0.0 );

    return( Runoff );

}


// private functions  ///////////////////////////////////////////////////////////
float Infiltration::NoInput( float Dt)     
{
  float DeepPerc = 0.0;

  if( Storage < FieldCapStorage ) {  // some deep perc and change in storage
    float diff_store = FieldCapStorage - Storage;
    float seepage = FinalRate*Dt;

    if( diff_store >= seepage){
      DeepPerc = seepage;
      Storage += seepage;
    }
    else {
      DeepPerc = diff_store;
      Storage = FieldCapStorage;
    }
  }

  return (DeepPerc);
}

float Infiltration::FillSoilToFC( float InFlux, float Dt )     
{

  float diff_store = Storage - FieldCapStorage;
  float InputVol = InFlux * Dt;
  float temp_dt = diff_store / InFlux;

  if( temp_dt >= Dt ) {  // not enough input to fill to FC
    Storage -= InputVol;
    temp_dt = Dt;
  }
  else {
    Storage = FieldCapStorage;
  }

  return( Dt - temp_dt );

}

float Infiltration::DrainSoilToTarget( float TargetStorage, float InFlux, float Dt, float *Perc )     
{

  float temp_dt = 0.0;
  float diff_store = TargetStorage -  Storage;
  float InputVol = InFlux * Dt;

  if( FinalRate == InFlux )
    temp_dt = Dt;
  else
    temp_dt = diff_store / (FinalRate - InFlux);
 
  if( temp_dt >= Dt ) {
    Storage -= ( InFlux - FinalRate ) * Dt;
    temp_dt = Dt;
  }
  else {
    Storage = TargetStorage;
  }

  *Perc = FinalRate * temp_dt;

  return( Dt - temp_dt );

}


// Interface function for C code --------------------------------------------
float CalcInfiltration( float store, float coef, float exp, float max_store, 
                        float fc_store, float final, float Dt, float InFluxVol, float *DeepPerc, 
                        float *RunoffVol  )
{
  float temp = 0, potential_store = 0;;
  // initial store = Moist[0], if greater than max_store, layer is super saturated

  // TEMP FIX - STILL NEED TO FIGURE CAUSE - MAY BE IN SURFACE ROUTING COMPONENT
  if(InFluxVol < 0){
    cout << "Negative InFluxVol supplied to CalcInfiltration: " << InFluxVol << endl;
    if( InFluxVol < -1e-10)exit(1);
    InFluxVol = 0.0;
  }

  if( store > max_store )
    temp = store - max_store;
  else
    potential_store = max_store - store; 
  
  Infiltration Surface(  potential_store, coef, exp, max_store, fc_store, final  );

  *RunoffVol = Surface.CalcRunoff( InFluxVol + temp, Dt, DeepPerc );

  // debug
  //float new_store = max_store - Surface.get_storage();
  //float error = new_store - store - InFluxVol + *RunoffVol + *DeepPerc;

  
  return ( max_store - Surface.get_storage() );
}



// test code ------------------------------------------------------------
#ifdef TEST

#include <iostream>
#include <fstream>

int main()
{
  int i = 0;
  float Storage = 0.0, DeepPerc = 0.0, RunoffVol = 0.0, OldStorage = 0.0, Error = 0.0;
  float RootDepth = 0.10, Porosity = 0.50, FieldCap = 0.20;
  float Coeff = 1.5, exp = 1.5, final = 0.002, Dt = 1.0, Rain = 0.006; //0.01;

  // open input file
  ifstream rain_file( "mica_rain.dat" );
  if( !rain_file.is_open() ) {
    cout << "could not open rain file " << endl;
    exit(1);
  }

  // open output file
  ofstream out_file( "output.test" );
  if( !out_file.is_open() ) {
    cout << "could not open output file " << endl;
    exit(1);
  }

  Infiltration Surface( 0.5*RootDepth, Coeff, exp, Porosity * RootDepth, FieldCap * RootDepth, final );

  for( i = 0; i < 5000; i++ ){
    
    //if( i > 25 ) Rain = 0.0;

    rain_file >> Rain;

    OldStorage = Surface.get_storage();
    RunoffVol = Surface.CalcRunoff( Rain*Dt, Dt, &DeepPerc );
    Storage = Surface.get_storage();
    Error = Rain * Dt + (Storage - OldStorage) - DeepPerc - RunoffVol;
        
    // print results
    //cout << "Storage " << Storage << "  Perc " << DeepPerc 
    //     << "  Runoff " << RunoffVol << "  Error " << Error << endl;

    out_file << Porosity*RootDepth - Storage << "  " << DeepPerc 
         << "  " << RunoffVol << "  " << Error << endl;

  }


  rain_file.close();
  out_file.close();

  return(0);


}
#endif


