#ifndef __FILE_SQUASH_H_SEEN__
#define __FILE_SQUASH_H_SEEN__

/*-----------------------------------------------------------------------------

Copyright (C) 2004, 2006.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program 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 2 of the License, or
(at your option) any later version.

This program 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 this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-----------------------------------------------------------------------------*/

#include "libscl.h"

namespace {

  const REAL pi = 3.14159265358979312e+00;
  const REAL a = pi/4.0;
  
  inline 
  REAL x0(REAL x, REAL c) 
  { 
    return x;
  }

  inline 
  REAL x1(REAL x, REAL& dydx, REAL c) 
  {
    dydx = 1.0; 
    return x;
  }

  inline 
  REAL x2(REAL x, REAL& dydx, REAL& dydx2, REAL c) 
  {
    dydx = 1.0; 
    dydx2 = 0.0; 
    return x;
  }

  inline 
  REAL arc0(REAL x, REAL c) 
  {
    if (x < -c) return std::atan(a*(x + c))/a - c;
    if (x >  c) return std::atan(a*(x - c))/a + c;
    return x;
  }
  
  inline 
  REAL arc1(REAL x, REAL& dydx, REAL c) 
  {
    if (x < -c) {
      REAL u = a*(x + c);
      dydx = 1.0/(1.0 + u*u);
      return std::atan(u)/a - c;
    }
    if (x >  c) {
      REAL u = a*(x - c);
      dydx = 1.0/(1.0 + u*u);
      return std::atan(u)/a + c;
    }
    dydx = 1.0;
    return x;
  }
  
  inline 
  REAL arc2(REAL x, REAL& dydx, REAL& dydx2, REAL c) 
  {
    if (x < -c) {
      REAL u = a*(x + c);
      dydx = 1.0/(1.0 + u*u);
      dydx2 = -2.0*u*a*dydx*dydx;
      return std::atan(u)/a - c;
    }
    if (x >  c) {
      REAL u = a*(x - c);
      dydx = 1.0/(1.0 + u*u);
      dydx2 = -2.0*u*a*dydx*dydx;
      return std::atan(u)/a + c;
    }
    dydx = 1.0;
    dydx2 = 0.0;
    return x;
  }
  
  inline
  REAL ave0(REAL x, REAL c) 
  {
    return 0.5*(arc0(x,c) + x0(x,c));
  }
  
  inline
  REAL ave1(REAL x, REAL& dydx, REAL c) 
  {  
    REAL dadx, dxdx;
    REAL y = 0.5*(arc1(x,dadx,c) + x1(x,dxdx,c));
    dydx = 0.5*(dadx + dxdx);
    return y;
  }
  
  inline
  REAL ave2(REAL x, REAL& dydx, REAL& dydx2, REAL c) 
  {  
    REAL dadx, dxdx;
    REAL dadx2, dxdx2;
    REAL y = 0.5*(arc2(x,dadx,dadx2,c) + x2(x,dxdx,dxdx2,c));
    dydx = 0.5*(dadx + dxdx);
    dydx2 = 0.5*(dadx2 + dxdx2);
    return y;
  }
  
  REAL gst0(REAL x, REAL c) 
  {
    REAL e = exp(x/c);
    REAL r = e/(1+e);
    REAL y = (4.0*c)*r - 2.0*c;
    return y;
  }
  
  REAL gst1(REAL x, REAL& dydx, REAL c) 
  {  
    REAL e = exp(x/c);
    REAL r = e/(1+e);
    REAL y = (4.0*c)*r - 2.0*c;
    REAL drdx = (r - r*r)/c;
    dydx = (4.0*c)*drdx;
    return y;
  }
  
  REAL gst2(REAL x, REAL& dydx, REAL& dydx2, REAL c) 
  {  
    REAL e = exp(x/c);
    REAL r = e/(1+e);
    REAL y = (4.0*c)*r - 2.0*c;
    REAL drdx = (r - r*r)/c;
    REAL drdx2 = (drdx - 2.0*r*drdx)/c;
    dydx = (4.0*c)*drdx;
    dydx2 = (4.0*c)*drdx2;
    return y;
  }
  
}

class squash {
private:
  REAL c;
  REAL (*f0)(REAL, REAL);
  REAL (*f1)(REAL, REAL&, REAL);
  REAL (*f2)(REAL, REAL&, REAL&, REAL);
public:
  squash() : c(), f0(&x0), f1(&x1),f2(&x2) { } 
  void set_identity() { f0 = &x0; f1 = &x1; f2 = &x2; }
  void set_logistic(REAL scale) {c=scale;; f0=&gst0; f1=&gst1; f2=&gst2;}
  void set_spline(REAL knot) {c=knot; f0=&ave0; f1=&ave1; f2=&ave2;}
  REAL operator()(REAL x) { return (*f0)(x,c); }
  REAL operator()(REAL x,REAL& dydx) { return (*f1)(x,dydx,c); }
  REAL operator()(REAL x,REAL& dydx,REAL& dydx2) {return (*f2)(x,dydx,dydx2,c);}
};

#endif
