// pntsrc_appx.hh
// Author: Mario J. Bencomo
// last modified 11/10/16

/**
 * \file pntsrc_appx.hh
 * For computing weights associated with singular source appx.
 */

#ifndef __PNTSRC_APPX_HH_
#define __PNTSRC_APPX_HH_

#include "MPS_includes.hh"

#define Q_SIZE 1+4+9+16+25+36

namespace TSOpt{
  
  static int Q_MAX  = 5;
  static float Q_MAT[Q_SIZE] = {
    // Q=0 -----------------------------------------------------------------------------
    1.0, 
    // Q=1 -----------------------------------------------------------------------------
    1.0, -1.0,
    0.0,  1.0,
    // Q=2 -----------------------------------------------------------------------------
    0.0, -0.5,  0.5,
    1.0,  0.0, -1.0,
    0.0,  0.5,  0.5,
    // Q=3 -----------------------------------------------------------------------------
    0.000000000000000,  -0.333333333333333,   0.500000000000000,  -0.166666666666667,
    1.000000000000000,  -0.500000000000000,  -1.000000000000000,   0.500000000000000,
    0.000000000000000,   1.000000000000000,   0.500000000000000,  -0.500000000000000,
    0.000000000000000,  -0.166666666666667,   0.000000000000000,   0.166666666666667,
    // Q=4 -----------------------------------------------------------------------------
    //row1
    0.000000000000000,   0.083333333333333,  -0.041666666666667,  -0.083333333333333,
    0.041666666666667,
    //row2
    0.000000000000000,  -0.666666666666667,   0.666666666666666,   0.166666666666667,
    -0.166666666666667,
    //row3
    1.000000000000000,   0.000000000000000,  -1.250000000000000,   0.000000000000000,
    0.250000000000000,
    //row4
    0.000000000000000,   0.666666666666667,   0.666666666666667,  -0.166666666666667,
    -0.166666666666667,
    //row5
    0.000000000000000,  -0.083333333333333,  -0.041666666666667,   0.083333333333333,
    0.041666666666667,
    // Q=5 -----------------------------------------------------------------------------
    //row1
    0.000000000000000,   0.050000000000000,  -0.041666666666667,  -0.041666666666667,
    0.041666666666667,  -0.008333333333333,
    //row2
    0.000000000000000,  -0.500000000000000,   0.666666666666666,  -0.041666666666667,
    -0.166666666666667,   0.041666666666667,
    //row3
    1.000000000000000,  -0.333333333333333,  -1.250000000000000,   0.416666666666667,
    0.250000000000000,  -0.083333333333333, 
    //row4
    0.000000000000000,   1.000000000000000,   0.666666666666667,  -0.583333333333333,
    -0.166666666666667,   0.083333333333333,  
    //row5
    0.000000000000000,  -0.250000000000000,  -0.041666666666667,   0.291666666666667,
    0.041666666666667,  -0.041666666666667, 
  //row6
    0.000000000000000,   0.033333333333333,   0.000000000000000,  -0.041666666666667,
    0.000000000000000,   0.008333333333333
  };
  
  /**
   * Truncated factorial: m!/(m-s)! = m*(m-1)*...(m-s+1)
   */
  int trunc_fact(int m,int s);
  
  /**
   * This function computes the ith (i=indx) weight in the pnt src approximation
   * by multiplying ith row of appropriate hard coded inverse matrix with some
   * appropriate RHS vector b.
   *
   * @param [in] s - derivative order
   * @param [in] p - approximation order
   * @param [in] indx -
   * @param [in] alpha - source grid residual
   * @param [in] h - spatial discretization
   */
  float pntsrc_weight( int s, int p, int indx, float alpha, float h);
  float pntsrc_weight( int sten_n, int s, int p, int indx, float alpha, float h);
}
#endif