package signal.lib;

import java.util.Vector;

/**
   Applies the Gegenbauer Reconstruction Procedure (GRP). Reconstruction parameters, m and L, are specified 
   locally.  The parameters are set independently in each smooth subinterval.  The intial parameters will be 
   the ones output from the postProcess() method of the gegenbauerReconstruction class.
*/


//   The values of m and L may than be changed in each interval via a dialog box in the GUI.
   


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

//           In
//  u      function passed to the method to be postprocessed
//  d      location of edges (with d[0]=ltBoundaryPt and d[nd+1]=rtBoundaryPt), created by edgeDetect or manually input 
//  mv     vector containing m in each subinterval which was created by gegenbauerReconstruction.postProcess()
//  Lv     vector containing L in each subinterval which was created by gegenbauerReconstruction.postProcess()

//           Out
//  uG     the postprocessed function that the method returns

 


   public static void postProcess(double u[], double d[], 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);   

// ---------------------------------------------------------------------------------   
// --- 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];
      Double Ld; Integer md;
      for (int j=0; j<(numDiscont-1); j++) {
         md = (Integer)mv.elementAt(j);  Ld = (Double)Lv.elementAt(j);  
         ma[j] = md.intValue();          La[j] = Ld.doubleValue();
      } // 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) ); }
              if (m==-1) uG[i] = u[i];
              else       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) ); }
              if (m==-1) uG[i] = u[i];
              else       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) ); }
                 if (m==-1) uG[i] = u[i];
                 else       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



} // gegenbauerReconstructB