package signal.lib;


import java.util.Vector;

   
/**
   Applies the Gegenbauer Reconstruction Procedure (GRP). Reconstruction parameters, m and L, are specified 
   globally through the parameters mK and LK.  After the postProcess method of this class is run on the function once,
   the postProcess() method in gegenbauerReconstructionB may be run and the reconstruction parameters may be set 
   separately in each interval.  The intial parameters will be the one output from the postProcess() method of this class.
   
   @see "On the Gibbs Phenomenon and its resolution, SIAM Review 39 (1998), p. 664-668"
*/

public final class gegenbauerReconstruction {
  private static utility SF;
  private static chebyshevPolynomial C;
  private static gegenbauerPolynomial G;

//       IN
//  u      function to be postprocessed
//  d      location of edges (with d[0]=ltBoundaryPt and d[nd+1]=rtBoundaryPt), created by edgeDetect or manually input
//  globally specify the reconstruction parameters
//  LK     Gegenbauer Exponent parameter in L = 0.5*LK*(d[j+1]-d[j])*N
//  mK     parameter in order of the Gegenbauer Polynomial m = 0.5*MK*(d[j+1]-d[j])*N

//       Out
//  uG     the postprocessed function that the method returns
//  mv     vector of m in subintervals (mv can then be passed to gegegenbauerReconstructionB as an initial choice of m in each subinterval)
//  Lv     vector of L in subintervals (Lv vector can then be passed to gegegenbauerReconstructionB as an initial choice of L in each subinterval)
   

   public static void postProcess(double u[], double d[], double LK, double MK,
                                  double uG[], Vector mv, Vector Lv) {
                 
    int        N = u.length - 1;                                         
    double  ak[] = new double[N+1];   // Chebyshev coefficients
    double   x[] = new double[N+1];
    int numDiscont = d.length;
    int sum=0;
    
     for ( int j=0; j<=N; j++) x[j] = C.xGLr(j,N);   // reverse order, i.e, 1.0...-1.0

// ---------------------------------------------------------------------------------   
// --- obtain Chebyshev coefficients for the function on the whole interval [-1,1] -
// ---------------------------------------------------------------------------------
 
    C.chebyCoef(u,ak,N);

// ---------------------------------------------------------------------------------
// ------  the reconstructed function  as  ug(i)= sum[ fh(L,m)*C(L,m,x) ] ----------
// ---------------------------------------------------------------------------------


      double L=0;  int m=0;
      int[] ma = new int[numDiscont];   double[] La = new double[numDiscont];
      
      for (int j=0; j<(numDiscont-1); j++) {
         m = (int)Math.ceil(0.5*MK*(d[j+1]-d[j])*N);   L = 0.5*LK*(d[j+1]-d[j])*N;
         mv.addElement(new Integer(m));                Lv.addElement(new Double(L));
         ma[j]=m;                                      La[j]=L;
      } // j

      double xi=0,K=0;
      
      int subinterval=0, oldSubInt=-1;
      for (int i=0; i<=N; i++) {
         uG[i]=0; xi=x[i];
         
         if (i==N) { // xi = -1.0
            subinterval=0; m = ma[0]; L = La[0];
            if (subinterval!=oldSubInt) { oldSubInt=subinterval; K = G.km(m,L)/( G.km(m+1,L)*G.h(m,L) ); }
            uG[i] = gh(ak,m,L,N,d[0],d[1],xi,K);
         } else if (i==0) { // xi = 1.0
             subinterval=numDiscont-2; m = ma[numDiscont-2]; L = La[numDiscont-2];
             if (subinterval!=oldSubInt) { oldSubInt=subinterval; K = G.km(m,L)/( G.km(m+1,L)*G.h(m,L) ); }
             uG[i] = gh(ak,m,L,N,d[numDiscont-2],d[numDiscont-1],xi,K);
         } else {
         
         
           for (int j=0; j<(numDiscont-1); j++) {
              if (xi>=d[j] && xi<d[j+1]) {
                subinterval=j; m = ma[j]; L = La[j];
                if (subinterval!=oldSubInt) { oldSubInt=subinterval; K = G.km(m,L)/( G.km(m+1,L)*G.h(m,L) ); }
                uG[i] = gh(ak,m,L,N,d[j],d[j+1],xi,K);
                break;    // get out of the j loop, the interval has been found
              } // if
           } // j
         
         } // else
         
      } // i
    } // gegenReconstruct
  
     

   private static double gh(double ak[], int m, double L, int N, double a, double b, double xi, double K) {
          
       int Nq = N + 2*(int)(L+1) + m;   
       double NI = 1.0/Nq;
       double sum=0, diff, gpL, g, gpm, gmmpL, gj, gpj, gp, xj, x;
       double eps = 0.5*(b-a), delta = 0.5*(b+a);
       
       x = G.tx(a,b,xi);
       
       for (int j=0; j<=Nq; j++) {
          xj = C.xGLr(j,Nq);
          diff = Math.abs(x-xj);
      
          if (diff<0.00000001) {    // if x == xj
          
                gpL  = NI*G.gegenPoly(m,L+1,xj);     g = Math.PI*K*G.gegenPoly(m,L,xj); 
                gpm  = NI*G.gegenPoly(m+1,L,xj);  gmmpL = Math.PI*K*G.gegenPoly(m-1,L+1,xj);
           
             sum += C.chebySeries(ak,eps*xj+delta)*Math.pow(1.0-xj*xj,L)*C.CJ(j,Nq)*2.0*L*(gpL*g-gpm*gmmpL);
      
          } else {
      
                gj  = NI*G.gegenPoly(m,L,xj);          gpj = NI*G.gegenPoly(m+1,L,xj); 
                 g  = Math.PI*K*G.gegenPoly(m,L,x);     gp = Math.PI*K*G.gegenPoly(m+1,L,x);
             
            sum += C.chebySeries(ak,eps*xj+delta)*Math.pow(1.0-xj*xj,L)*C.CJ(j,Nq)*( (gpj*g-gj*gp)/(xj-x) );
      
         } // if
      } // j
      
     return  sum;

  } // gh


} // gegenbauerReconstruction