/* 
* 
* Copyright (C) 2013 UMPC
* 
* This file is part of CLaRiNet.
*  
* CLaRiNet is free software: you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* CLaRiNet is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with CLaRiNet.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <stddef.h>
#include "GeoMorph.h"
#include "RasterCalc.h"
#include "Erreur.h"
#include "esriheader.h"
#include "FloatRaster.h"
#include "Topology.h"
#include "CellTree.h"
#include "Hypsometry.h"
#include <math.h>

using namespace std;

int SyntheticHypso ( const int HandleTree, const int HandleRaster , float smax , float alpha , float A0 , float zbase )
{
    CellTree * T;
    FloatRaster * R;

    if( !(T = GetCellTree(HandleTree) ) )
    {
        SetErreur("SyntheticHypso : tree # %d does not exist\n",HandleTree);
	return(-1);
    }

    if( !(R = GetFloatRaster(HandleRaster) ) )
    {
        SetErreur("SyntheticHypso : raster # %d does not exist\n",HandleRaster);
	return(-1);
    }

    if ( SyntheticHypso_ ( T , & (T->RootCell) , & zbase , & smax , & alpha , & A0 , R ) )
    {
        SetErreur("SyntheticHypso : recursive call failed\n");
	return(-1);
    }

    return(0);
}

int SyntheticHypso_ ( CellTree * T , int * _recDown , float * _zDown , float * _smax , float * _alpha , float * _A0 , FloatRaster * R )
{
    float zcell , znext , dL ,s;
    int d2 , recnext ;

    // On calcule l'élévation théorique par rapport à la cellule aval

    (void)  d2rec( & (R->Header) , _recDown , & (T->RootCell) , & d2 );
    dL = sqrt((float)d2);
    dL = ((R->Header).cellsize)*dL;

    s = ((float)(T->Size))/(*_A0);
    s = powf(s,(*_alpha));
    s = (*_smax)*s/(1+s);       // s = smax * (A/A0)^alpha / ( 1 + (A/A0)^alpha )

    zcell = (*_zDown) + dL*s; // z = z0 + s *dL

    (void) PutRasterValue1 ( R , & (T->RootCell) , & zcell );

    // Appels récursifs

    if ( T->FirstChild )
    {
        znext = zcell ;
        recnext = (T->RootCell);

        if ( SyntheticHypso_ ( T->FirstChild , & recnext , & znext , _smax , _alpha , _A0 , R ) )
        {
            SetErreur("SyntheticHypso_ : recursive call failed\n");
            return(-1);
        }
    }

    if ( T->NextSibling )
    {
        znext = (*_zDown) ;
        recnext = (*_recDown);

        if ( SyntheticHypso_ ( T->NextSibling , & recnext , & znext , _smax , _alpha , _A0 , R ) )
        {
            SetErreur("SyntheticHypso_ : recursive call failed\n");
            return(-1);
        }
    }

    return(0);

}
