#include "channels.h"
#include "rungekutta.h"

extern "C" {
float route_segment_powerlaw(float inFlux, float dt, float startDepth, float width, float length, float slope, 
                    float roughness, float *storage); 
} 


float ChannelWater::ds_inFlux;

// Abastract Base class for Channel  ------------------------------------------
ChannelWater::ChannelWater( float startDepth, float length, float slope, float roughness ) {
  channelStateWater( startDepth, length, slope, roughness );
}


void ChannelWater::channelStateWater( float startDepth, float length, float slope, float roughness ) {
  assert( startDepth >= 0 && length >= 0 );
  d_waterDepth = startDepth;
  d_length = length;
  d_slope = sqrt( slope );
  d_roughness = roughness;
}

float ChannelWater::calcDerivative( float dummy, float depth ){
  return (ds_inFlux - (d_slope/d_roughness)*crossSecArea(depth)*pow(hydRadius(depth), 2./3.) )
         / surfaceArea(depth);
}

float ChannelWater::solveBalance( float inVol, float startDepth, float endDepth ) {
    assert( inVol >= 0 );

    float outVol = inVol - ( storageVol(endDepth) - storageVol(startDepth) );
    if( outVol > 0 ){
      d_waterDepth = endDepth;
    }
    else{
      outVol = storageVol(startDepth) + inVol;
      d_waterDepth = 0;
    }

    return outVol;
}


// Final RectangularChannel Class -----------------------------------------------------------
RectangularChannelWater::RectangularChannelWater( float startDepth, float width, float length,
                                                 float slope, float roughness )
                                                 : ChannelWater( startDepth, length, slope, roughness ){
  d_bottomWidth = width;
}

void RectangularChannelWater::initialStateWater( float startDepth, float width, float length,
                                                 float slope, float roughness) {
  assert( width >= 0 );
  channelStateWater( startDepth, length, slope, roughness );
  d_bottomWidth = width;
}

inline float RectangularChannelWater::surfaceArea( float depth ){ 
  return d_bottomWidth * d_length; 
}

inline float RectangularChannelWater::crossSecArea( float depth ){ 
  return d_bottomWidth * depth; 
}

inline float RectangularChannelWater::hydRadius(float depth ){ 
  return (d_bottomWidth * depth) / (d_bottomWidth + 2*depth); 
}

inline float RectangularChannelWater::storageVol(float depth ){ 
  return crossSecArea( depth ) * d_length; 
}

#ifdef LINK_DHSVM


// Interface for C code  -----------------------------------------------------------------
float route_segment_powerlaw(float inFlux, float dt, float startDepth, float width, float length, float slope, 
                    float roughness, float *storage) 
{

  RectangularChannelWater RectChannel( startDepth, width, length, slope, roughness );


  RectChannel.getInflux( inFlux );
  float newDepth = RungeKutta( RectChannel, 0.25*dt, 0, dt, startDepth, 1e-5 );

  float outflow = RectChannel.solveBalance( inFlux *dt, startDepth, newDepth );

  *storage = RectChannel.storageVol( RectChannel.waterDepth() );

  return(outflow);
}

#endif





