// PowerLaw.cpp

// Solves general ODE: dS/dt = I - c(S^n)

#include <iostream>
#include <cstring>
#include <fstream>
#include <cassert>
#include <cmath>
#include "powerlaw.h"

#define TEST
#undef TEST

// public functions
PowerLaw::PowerLaw()
  {
    Storage = 0.0;
    Coeff = 0.0;
    Exponent = 1.0;
  }

PowerLaw::PowerLaw( float store, float coef, float exp )
  {
    assert( store >= 0.0 );
    Storage = store;
    Coeff = coef;
    Exponent = exp;
  }

void PowerLaw::get_parameters( float store, float coef, float exp )     // make virtual 
{
    assert( store >= 0.0 );
    Storage = store;
    Coeff = coef;
    Exponent = exp;
}



float PowerLaw::output( float InputVol, float Dt )     // make virtual 
{
  assert( InputVol >= 0 );
  float NewStorage = 0.0;
  float InFlux = InputVol/Dt;

  if( Exponent == 1.0 ) // linear - analytical solution
  {
    NewStorage = CalcAnalytical( InFlux, Coeff, Dt );
  }
  else  // non-linear powerlaw
  {
    NewStorage = ODEsolver( InputVol/Dt, Storage, Dt );
  }

  if( NewStorage < 0.0 ) NewStorage = 0.0;
  float Output = InputVol - ( NewStorage - Storage );
  Storage = NewStorage;

  return( Output );

}


// private/protected functions ------------------------------------------

// calculate the derivative (unique for power law storage relationship ) //////
float PowerLaw::CalcDerivative( float Influx, float NewStore )  
{
  if (NewStore <= 0.0) {return (Influx);} // SRW 4-20-01.  Linux doesn't like pow() when argument 1 is 0.0!
  else {return( Influx - Coeff * pow( (double) NewStore, (double) Exponent ) );}
}


// private base functions /////////////////////////////////////////////////////

// calculate analytical solution if Exponent equals 1.0  //////////////////////
float PowerLaw::CalcAnalytical( float InFlux, float k, float Dt )
{
  return( (InFlux/k)+( Storage - (InFlux/k) )*exp(-k*Dt) );
}

// calculate new storage ( general for dS = I - Q(S) )
// using RungeKutta
float PowerLaw::ODEsolver( float InFlux, float Storage, float Dt )
{
  float NewStorage, k1, k2, k3, k4;

  k1 = Dt * CalcDerivative( InFlux, Storage );
  k2 = Dt * CalcDerivative( InFlux, Storage + 0.5*k1 );
  k3 = Dt * CalcDerivative( InFlux, Storage + 0.5*k2 );
  k4 = Dt * CalcDerivative( InFlux, Storage + k3 );
  NewStorage = Storage + (1.0/6.0) * ( k1 + 2*k2 + 2*k3 + k4 );

  if( NewStorage < 0.0 ) NewStorage = 0.0;  // SRW 4-17-01 

  return ( NewStorage );
}

// test code ------------------------------------------------------------
#ifdef TEST
int main()
{
  float Ks = 1.0;                 // m/h
  float Slope = 0.3;              // m/m
  float Exponent = 1.0;
  float InFlux = 0.05;            // m/h
  float Dt = 1.0;

  float Outflow = 0.0, Storage = 0.0;
  float OutFlowA = 0.0, StoreA = 0.0, NewStoreA = 0.0, k = 0.0;
  PowerLaw Sat( Storage, Ks * Slope, Exponent );      // create Sat

  k = Ks * Slope;


  // loop through
  for( int i = 0; i < 10; i++ )
  {
    Outflow = Sat.output( InFlux * Dt, Dt );   // RungeKutta
    Storage = Sat.get_storage();

    // analyitical w/ exponent = 1 (linear storage resv. )
    NewStoreA = (InFlux/k)+(StoreA - (InFlux/k))*exp(-k*Dt);
    OutFlowA = InFlux*Dt - ( NewStoreA - StoreA );

    // print results
    cout << "Out " << Outflow << "  OutA " << OutFlowA 
         << "  Store " << Storage << "  StoreA " << NewStoreA << endl;

    StoreA = NewStoreA;
  }

  return( 0 );

}

#endif


