#include "scltypes.h"
#include "sclfuncs.h"
#include "realmat.h"

using namespace std;
using namespace scl;


#define MAXK 100

int main() 
{

  REAL h; 
  REAL x;  
  INTEGER K=10;  
  REAL c[MAXK+1];  
  REAL dc[MAXK+1];
  REAL s[3];

  INTEGER i;
  REAL sum, tmp;

  ofstream fout("tsthermi.plt");  //plot output for scatter
  if (!fout) {
    cerr << "Error, tsthermi, cannot open file tsthermi.plt\n" ; 
    exit(1);
  }

  for (sum=0.0, i=0; i<=K; i++) {
    c[i]=1.0/pow(REAL(i+1),4); 
    sum += pow(c[i],2);
  }
  tmp = 1.0/sqrt(sum);
  for (i=0; i<=K; i++) {
    c[i] *= tmp; 
  }

  for (sum=0.0, i=0; i<=K; i++) {
    sum += pow(c[i],2);
  }

  cout << starbox("/Sum of squares of c//") << sum << '\n';

  cout << starbox("/c//") ;
  dgmpnt(cout,c,K+1,1) << starbox("/Dgmpnt returns a reference to cout.//");

  REAL xl = -4.0; 
  REAL xu =  4.0;

  REAL xinc;

  xinc = (xu-xl)/(LINESIZE-1.0);

  h = hermite(0.0, 0, c, dc, s);

  if (fabs( (double)(h/(*c) -(*dc)) ) > 1.e-10)
     cout << starbox("/Grevious error///") <<  h/(*c) << "  " << (*dc) << '\n';

  for (x = xl; x <= xu; x += xinc) {
    h = hermite(x, K, c, dc, s);
    #if defined GNU_GPP_COMPILER || defined HP_CC_COMPILER
      fout << "  " << h << "  " << x << '\n';
    #else
      fout << setw(25) << h << setw(25)  << x << '\n';
    #endif
  }

  xinc = 1.e-2;

  for (sum = 0.0, x = xl; x <= xu; x += xinc) {
    h = hermite(x, K, c, dc, s);
    sum += h*h*xinc;
  }

  cout << starbox("/integral of SNP //") << sum << '\n';

  cout << starbox("/derivative check//") << '\n';

  #if defined GNU_GPP_COMPILER || defined HP_CC_COMPILER
    cout << "   " << "i" << "     " << "dc[i]"
      << "     " << "(hplus-h)/eps" << '\n';
  #else
    cout << setw(30) << "i" << setw(15) << "dc[i]"
      << setw(15) << "(hplus-h)/eps" << '\n';
  #endif

  REAL eps = 1.e-8;
  REAL hplus, csave;

  x=0.0;

  for (i=0; i<=K; i++) {
    csave = c[i];
    c[i] += eps;
    hplus = hermite(x, K, c, dc, s);
    c[i] = csave;
    h = hermite(x, K, c, dc, s);
    #if defined GNU_GPP_COMPILER || defined HP_CC_COMPILER
      cout << "   " << i << "     " << dc[i]
        << "     " << (hplus-s[0])/eps << '\n';
    #else
    cout << setw(30) << i << setw(15) << dc[i]
      << setw(15) << (hplus-s[0])/eps << '\n';
    #endif
  }

  cout << '\n';

  REAL xsave = x;  REAL splus[3];
  x += eps;
  hplus = hermite(x, K, c, dc, splus);
  x = xsave;
  h = hermite(x, K, c, dc, s);

  #if defined GNU_GPP_COMPILER || defined HP_CC_COMPILER
    cout << "       " << "s[1]"  << "       " << "(splus[0]-s[0])/eps" << '\n';
  #else
    cout << setw(25) << "s[1]"  << setw(25) << "(splus[0]-s[0])/eps" << '\n';
  #endif

  #if defined GNU_GPP_COMPILER || defined HP_CC_COMPILER
    cout << "       " << s[1] << "       " << (splus[0]-s[0])/eps << '\n';
  #else
   cout  << setw(25) << s[1] << setw(25) << (splus[0]-s[0])/eps << '\n';
  #endif

  cout << '\n';

  #if defined GNU_GPP_COMPILER || defined HP_CC_COMPILER
    cout << "     " << "s[2]" << "       " << "(splus[1]-s[1])/eps" << '\n';
  #else
    cout << setw(25) << "s[2]" << setw(25) << "(splus[1]-s[1])/eps" << '\n';
  #endif

  #if defined GNU_GPP_COMPILER || defined HP_CC_COMPILER
    cout << "       " << s[2] << "       " << (splus[1]-s[1])/eps << '\n';
  #else
  cout  << setw(25) << s[2] << setw(25) << (splus[1]-s[1])/eps << '\n';
  #endif


  fout.close();
  exit(0);
}


