/*====================================================================================
! Purpose(s): 
!
! Contain(s): CalcFrictionVelocity
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!====================================================================================*/

#include "CCalcFricVel.h"
#include "CErrors.h"
#include "math.h"
#include "CMathMisc.h"
#include "stdio.h"

double CalcFrictionVelocity(double AE, double AVt, double ADist, double ANu, double AKappa, bool AIsRoughSurface, double AKs) 
{
/*!=======================================================================
! Purpose(s): Calculates friction velocity. If wall is a smooth one assumes a log profile. If it
!             is a rough surface uses a different formulation.
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!=======================================================================*/ 
  double A, X, Xo, LogXo;
  double Value=0.0;
  int n;


  const double MaxErr=1.0e-3L;
  const int MaxIt=10;

  if (fabs(AVt)>1e-6)
  // if velocity is big enough...
  {
    if (AIsRoughSurface)
    // If it is a rough surface...
    {
      Value=AKappa*AVt/log(30.0*ADist/AKs);
    }
    else
    // If it is not a rough surface...
    {

      A=AE*fabs(AVt)*ADist/ANu;
      
      X=1.0;
      n=0;

      if (A>1.0)
      // For velocities large enough...
      {
        A=log(A);

        // We want to solve the equation: x·(ln(x)+A)-Akappa = 0.
        // For doing this we use an iterative Newton Raphson method.

        // function x·(ln(x)+A)-Akappa has only one minimum at x0=exp(-A-1). This minimun is negative,
        // and the function at the right of x0 is growing monotone, being the solution of x·(ln(x)+A)-Akappa=0
        // always at the right of this minimum.
        // As an inital guess we use 2*x0.
        X=2.0*exp(-A-1.0);

        Xo=0.0000001e0L;
        while (fabs((X-Xo)/Xo)>MaxErr && n<=MaxIt && X>0.0)
        {
          Xo=X;

          LogXo=log(Xo);

          // Usual Newton-Raphson's method.
          X = Xo + (AKappa-Xo*(A+LogXo))/(1.0+A+LogXo);
          
          // Extended Newton-Raphson using parabolas instead of lines.
          //X = DSqrt(Xo)*DSqrt(Xo*(LogXo*LogXo) + 2.0*Xo*A*LogXo + Xo*(A*A + 1.0) + 2.0*AKappa) - Xo*LogXo - Xo*A;

          n=n+1;
        } //while

        if (n>=MaxIt) 
        // If the number of iterations has reached the maximum...
        {

          printf("WARNING: too many iter", __FILE__, __LINE__);//dsqrt(ANu/ADist*dabs(AVt)), X*dabs(AVt)
        } //if

        Value = max2(X*fabs(AVt), sqrt(ANu/ADist*fabs(AVt)));
      
      }
      else
      // For really small velocities (that yields to negative logarithms)...
      {
        Value = sqrt(ANu/ADist*fabs(AVt));

      } //if
      
      if (X<=0 || n==MaxIt || Value<1e-6) 
      {
        //Error("****ERROR in CalcFrivCel", __FILE__, __LINE__); //x, n, Value
      }

      // Fric vel has the same sign as AVt.
      Value = copysign(Value, AVt);

    } //if

  }
  else
  // if velocity is small, frictional velocity is assumed to be 0.
  {
    Value = 0.0;
  } //if

  return Value;

};