#include <except.hh>
#include <su.h>
#include <segy.h>
#include <header.h>
#include <parser.h>

#include "MPS_includes.hh"
#include "pntsrc_appx.hh"

#ifndef TOL
#define TOL 1.e-7
#endif

const char * sdoc[] = { 
  "Usage: test_pntsrc_appx.x a_ord= d_ord= h= a=",
NULL};  

using TSOpt::pntsrc_weight;
using TSOpt::trunc_fact;
using RVL::valparse;

int xargc;
char **xargv;

//-----------------------------------------------------------------------------------
int main(int argc, char ** argv) {
//-----------------------------------------------------------------------------------
  try {
 
    xargc=argc; xargv=argv;
    requestdoc(1);

    PARARRAY * par = ps_new();
    if ( ps_createargs(par, argc - 1, argv + 1) ) {
      printf("Error parsing input data. ABORT.\n");
      exit(1);
    }
    ps_printall(*par,stderr);

    //parameters    
    float h = valparse<float>(*par,"h");
    float a = valparse<float>(*par,"a"); 
    int a_ord = valparse<int>(*par,"a_ord");
    int d_ord = valparse<int>(*par,"d_ord");
   
    //grid point info of approx src point
    int st_N = a_ord + d_ord;

    if(st_N<=0){
      RVLException e;
      e << "ERROR from test_pntsrc_appx.x: stencil size of pnt appx is <=0!";
      e << "st_N = " << st_N <<"\n";
      throw e;
    }
      
    vector<float> stencil;
    int j = (int)((st_N-1)/2);
    cerr << "j=" << j <<"\n";

    for( int i=0; i<st_N; i++){      
      //computing weights
      stencil.push_back(pntsrc_weight(d_ord,a_ord,i,a,h));     
      cerr << "stencil["<< i <<"]="<< stencil[i] <<"\n";
    }

    int order=0;
    //testing weights satisfy moment conditions
    for( int m=0; m<10; m++){
      
      //computing mth moment
      float mom = 0.0;
      for( int i=0; i<st_N; i++){
	 mom += pow(h*(i-j-a),m) * stencil[i];
      }
      mom *= h;
      cerr << "m=" << m <<" moment = " << mom <<"\n";

      float diff=1;
      //checking moment condition
      if(m==d_ord){
	mom *= pow(-1,d_ord)/trunc_fact(d_ord,d_ord);
	diff = fabs(mom-1.0);
	cerr << "diff=" << diff <<"\n";
      }
      else diff = fabs(mom);

      if(diff>TOL && order==0) order=m; 
    }
    
    order -= d_ord;
    cerr << "approximation order = " << order << "\n";

  }
  catch (RVLException & e) {
    e.write(cerr);
    exit(1);
  }
}