package signal.lib;


import java.util.Vector;

/**
 Finds edge locations in functions.  Only searches the interval [A,B].
*/

 public final class edgeDetectAB {
    private static chebyshevPolynomial CM;
    
    
// Same functionality as edgeDetect except only looks for edges on the
// subinterval [A,B]. The findEdges method is used to find edges in the
// derivative of the data in an interval [A,B] when the data has discontinuities
// in the first derivative (rarefaction waves).

//       IN
//  u          the function to be postprocessed, this data may be known on the 
//             CGL grid on [-1,1] or on a coord. trans. of this grid
//  xm         the grid that the function is know on
//  JCRIT      critical threshold
//  NE         if an edge is located at x[j], then only one edge may be detected 
//             the interval (x[j]-x[j-NE],x[j]+x[j+NE])
//  Q          enhancement exponent

//       OUT
//  d      vector containing the location of the edges in order -1,..,1
//         with d[0]=-1.0 and d[M]=1.0; thus, there will be (M-1) edges found
//  uE     edge series
//  uNLE   enhanced edge series


 public static void 
findEdges(double[] u, double JCRIT, int Q, double[] uE, double[] uNLE, Vector d, int NE, double A, double B, double[] xm) {
      int        N = u.length - 1;
      double    Nd = (double)N;
      double  aN[] =new double[N+1];              // Chebyshev coefficients
      double   x[] =new double[N+1];

       for ( int j=0; j<=N; j++) x[j]= Math.cos(Math.PI*j/Nd);   // order 1,..,-1

       
// --- find chebyshev coefficients ------------------------------

     double dx=(1.0/N), dxP=dx*Math.PI;
     double sum2=0;

     for (int k=0; k<=N; k++ ) {
       sum2=0;
       for (int n=0; n<=N; n++ )   sum2 += CM.CJ(n,N)*u[n]*Math.cos(k*n*dxP);
       aN[k]=2.0*CM.CJ(k,N)*dx*sum2;
     } // k 

// --------------------------------------------------------------
       

       for (int j=1; j<N; j++) {
         for (int k=1; k<=N; k++)   uE[j] +=  aN[k]*CM.Tp(x[j],k);
         uE[j] *= (Math.sqrt(2.0)*Math.PI*Math.sqrt(1.0-x[j]*x[j])/Nd);
       } // j

       uE[0]=0;
       uE[N]=0;
                                                                      // nonlinear enhancement
       double temp=0;
       for (int j=0; j<=N; j++) {
        temp = Math.abs( Math.pow(N,0.5*Q)*Math.pow(uE[j],Q) );
         if (temp<JCRIT) uNLE[j] = 0;                             // the direction of the jump is not important
                                                                  //  so take abs val
         else            uNLE[j] = Math.abs(uE[j]);
        } // j
        
     Vector tempV = new Vector();
     
//  ---  limit one edged per neighborhood of size 2*NE ---

     boolean addEdgeL=false, addEdgeR=false;
                                                                      // find edges
     for (int j=1; j<N; j++) {
     
      if (xm[j]>A && xm[j]<B) {
     
       if ( (j-NE)>=0 && (j+NE)<=N ) {
       
         if (uNLE[j]>0) {                                              // check surronding points
           for (int i=1; i<=NE; i++) {
              if (uNLE[j-i]>0 && (uNLE[j-i]<uNLE[j])) uNLE[j-i]=0;
              if (uNLE[j+i]>0 && (uNLE[j+i]<uNLE[j])) uNLE[j+i]=0;
           }
         }
       }  
        else {  // it is within NE grid points of a boundary
          uNLE[j]=0;
       }  // else
       
         } // if in (Ae,Be)
         else uNLE[j]=0; // x not in (Ae,Be)
 
       } // j
       
        for (int j=1; j<N; j++) if (uNLE[j]>0) tempV.addElement(new Double(x[j]));
     
     
// ----------------------
     
     // put edge location in correct order, -1,...,1
     int sizeTempV = tempV.size()-1;
     for (int j=0; j<=sizeTempV; j++) d.addElement(tempV.elementAt(sizeTempV-j));


 } // findEdges

} // edgeDetectAB