/* Last edited on 2007-01-08 02:05:21 by stolfi */ 

/*
 * Functions to compute with sparse decompositions of(E and H).
 *
 * The total number of cells n is specified indirectly, using the number of
 * detail levels, ell, and the number of cells at the coarsest approximation
 * level, n0.  It follows that n = n0*2^ell.
 * The variables n, n0, and ell used below always have this meaning.
 *
 * It is tacitly assumed that n0 is at least equal to the number coefficients
 * of the interpolating polynomials (that is, n0 >= two for linear
 * interpolation, and n0 >= 4 for cubic interpolation).
 */

  /* Sonia: validade do indice deve ser testada ANTES de usar o indice. [stolfi] */
  /* Os operadores && e || nćo sćo comutativos. [stolfi] */


/* We need to define _GNU_SOURCE because we use the constant NAN. */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdarg.h>

#include <int2d.h>


/*
 * LOCAL function
 *
 * Returns the decomposition level to which the k-th and m-th data sample belongs.
 */
void find_level(int *inf, int k, int m, int ell)
{
  int level=ell; /* num detail levels */
  inf[0]=0;
  inf[1]=0;
  inf[2]=0;
  while (level)
    {
      if((k & 1)||(m & 1))
        { 
          inf[0]=level;
          inf[1]=(k & 1);
          inf[2]=(m & 1);
          return;
        }
      else
        {         
          k >>= 1;
          m >>= 1;
          level--;
        } 
    }
  inf[0]=0;
  inf[1]=0;
  inf[2]=0;

}



/*=======================================================================*/
/*=======================================================================*/

/* LOCAL function
 *
 * This function does the basic interpolation of a function with
 * symmetric boundary conditions on the left and right sides, and anti
 * symmetric boundary contitions on the front and back sides. which is
 * the typical boundary conditions for Ex and Dy(Hz). It computes the
 * value u[k], given the sparse representation stored in u[0..n]. It
 * is tacitly assumed that u[k] is not in the sparse decomposition,
 * that is, u[k] is NAN. The number of points used to interpolate is
 * interp_points, and the spacing at the next coarsest level is spc.
 */
double interpEx(double** u, int ell, int n0, int interp_points, int k, int m)
{
  int n = (n0) * (1 << ell);
  double v[interp_points];
  double h[interp_points];
  int ix[interp_points];
  int iy[interp_points];
  double r;
  int spc;
  int inf[3];
  int j;
  int level;
  int fx;
  int fy;
  find_level(inf, k, m, ell);
  level=inf[0];
  fx=inf[1];
  fy=inf[2];
  spc= 1 << (ell - level + 1);


  /*memset(v, 0, sizeof(v));
    memset(h, 0, sizeof(h));
    memset(ix, 0, sizeof(ix));
    memset(iy, 0, sizeof(iy));*/
  /*
   * The (interior) points to interpolate from are the following:
   *
   *
   * if interp_points == 4:    (need to redefine it close to the boundaries)
   *
   *    ix[0] = k - 3*spc/2;
   *    ix[1] = k - spc/2;
   *    ix[2] = k + spc/2;
   *    ix[3] = k + 3*spc/2;
   *
   * The following loop handles the general case (interp_points even and
   * greater than two).
   */

        
  if ((fx==1) && (fy==0))     /* interpolaĆ§Ć£o segundo o eixo dos XX */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }

      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          //ix[1] = ix[1];
          //ix[2] = ix[2];
          //ix[3] = ix[3];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          //ix[0] = ix[0];
          //ix[1] = ix[1];
          //ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de E\n");
              exit(-1);
            }
          if (isnan(u[ix[j]][m]))
            {
              v[j] = interpEx(u, ell, n0, interp_points, ix[j], m);
            }
          else
            {
              v[j] = u[ix[j]][m];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j]; /*faz as contas*/
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx */


  else if ((fx==0) && (fy==1))     /* interpolaĆ§Ć£o segundo o eixo dos YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          iy[interp_points/2-j] = m - (2*j-1)*spc/2;
          iy[interp_points/2-1+j] = m + (2*j-1)*spc/2;
        }

      if (iy[0] < 0) /* left boundary */
        {
          iy[0] = -iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[3];
          h[0] =  1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (iy[3] > n) /* right boundary */
        {
          iy[0] = iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (iy[j] < 0 || iy[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, iy[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
          if (isnan(u[k][iy[j]]))
            {
              v[j] = interpEx(u, ell, n0, interp_points, k, iy[j]);                 
            }
          else
            {
              v[j] = u[k][iy[j]];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos yy*/

  else if ((fx==1) && (fy==1))   /* interpolaĆ§Ć£o segundo o eixo dos XX e YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }

      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[3];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          ix[0] = ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
                
          v[j] = interpEx(u, ell, n0, interp_points, ix[j], m); 
                 
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
           
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx e dos yy */
  else
    { printf("Estou na malha mais grossa. o valor de fx e %d e de fy e %d\n", fx, fy);
      exit(-1);
    }
}
/*=======================================================================*/
/*
 * PUBLIC function
 *
 * Converts a data vector u of type Ex  to a sparse representation.
 * It sets a data value u[i][m]  to NAN if the absolute value of the
 * corresponding wavelet coefficient  is smaller than epsilon.
 */
void spr_data_to_sparseEx(double **u, int ell, int n0, int interp_points, double eps)
{
  int n = (n0) * (1<<ell);
  int spc = 1;
  int i,j,m;


  /* run through the data, from the finest to the coarsest detail levels */
  for (j=ell; j >= 1; j--)
    {
      spc <<= 1; /* the separation between points doubles with each level */
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpEx(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      //printf("a fiferenca e %g\n", u[i][m] - interp(u, ell, n0, interp_points, i ,m, inf));
                      //scanf("%d",&p); 
                      u[i][m] = NAN;
                           
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas linhas*/


      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpEx(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpEx(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/


    }/*termina o ciclo em j*/
}

/*=======================================================================*/

void spr_sparse_to_dataEx(double **u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1<<ell);
  int spc =  1 << (ell);
  int i,j,m;
        

  /* run through the data, from the finest to the coarsest detail levels */
  for (j=1; j <= ell; j++)
    {
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (isnan(u[i][m]))
                u[i][m] = interpEx(u, ell, n0, interp_points, i ,m);
              //printf("OK\n");
            }
        }/*terminou a interpolaĆ§Ć£o nas linhas*/
                  
      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interpEx(u, ell, n0, interp_points, i ,m);
                         
            }                 
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interpEx(u, ell, n0, interp_points, i ,m);
                           
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/
      spc >>= 1; /* the separation between points doubles with each level */

    }/*termina o ciclo em j*/
}

/*=======================================================================*/
/* Recall that Ex is antisymmetric in y*/

void derivYEx(double** u,int ell, int n0, int interp_points, double dx, double **du,  double** v)
{

  int n = (n0) * (1<<ell);
  int i, m;

  for (i=0; i <=n; i++)
    for (m=0; m <=n; m++)   
      { if (isnan(v[i][m]))
          du[i][m] = NAN;
        else
          du[i][m]=spr_derivYEx(u, ell, n0, interp_points, dx, i, m);
      }

}

/*==========================================================================*/
double spr_derivYEx(double** u, int ell, int n0, int interp_points, double dx, int k, int m)
{
  int n = (n0) * (1 << ell);
  int ix[interp_points];        /* the positions of the samples needed to compute the derivative */
  double h[interp_points];      /* the weights needed to compute the derivative */
  double v[interp_points];      /* the function values needed to compute the derivative */
  int i;
  double du;

        
  if (m==0)
    {
      ix[0] = 2;
      ix[1] = 1;
      ix[2] = 1;
      ix[3] = 2;
      h[0] = -1./12.;
      h[1] = 2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (m==1)
    {
      ix[0] = 1;
      ix[1] = 0;
      ix[2] = 2;
      ix[3] = 3;
      h[0] = -1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (m==n-1)
    {
      ix[0] = n-3;
      ix[1] = n-2;
      ix[2] = n;
      ix[3] = n-1;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = 1./12.;
    }
  else if (m==n)
    {
      ix[0] = n-2;
      ix[1] = n-1;
      ix[2] = n-1;
      ix[3] = n-2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =-2./3.;
      h[3] = 1./12.;
    }
  else
    {
      ix[0] = m-2;
      ix[1] = m-1;
      ix[2] = m+1;
      ix[3] = m+2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] =-1./12.;
    }

 

  for (i=0; i <interp_points; i++)
    {
      if (ix[i] < 0 || ix[i] > n)
        {
          fprintf(stderr, "ERROR: spr_deriv(): ix[%d]=%d\n", i, ix[i]);
          printf("Estou na derivada de  Ex\n");
          exit(-1);
        }
      v[i] = u[k][ix[i]];
      if (isnan(v[i]))
        {
          v[i] = interpEx(u, ell, n0, interp_points, k, ix[i]);         /* interpolate u[] at ix[i] */  
        }
    }
  /* compute the derivative */
  du = 0;
  for (i=0; i < interp_points; i++) du += v[i] * h[i];
  du /= (1*dx);
  return du;
}

/*=======================================================================*/

//extern void** dvector(int,int);

void spr_refineEx(double** u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1 << ell);
  int i, m, i1 ,i2, i11, i21, m1, m2, m11, m21;
  int inf[3];
  double esp[n+1][n+1];
        
                

  for (i=0; i <=n; i++)
    for (m=0; m <=n; m++)
      esp[i][m]= u[i][m];

  for (m=0; m <= n;m++ )
    for (i=0; i <= n;i++ )
      {
        find_level(inf, i, m, ell);
        if ((!isnan(esp[i][m])) && (inf[0]>0))
          {
            i1= i-(1<<(ell-inf[0]));
            i2= i+(1<<(ell-inf[0]));                
            m1= m-(1<<(ell-inf[0]));
            m2= m+(1<<(ell-inf[0]));
            //printf("o valor de i1 e %d de i2 e %d de m1 e %d e de m2 e %d\n", i1, i2, m1, m2);
                 
            if ((i1 >= 0) && (isnan(esp[i1][m])))
              {
                u[i1][m] = interpEx(u, ell, n0, interp_points, i1, m);
              }
            if((i1 >= 0)  && (m2 <= n) && (isnan(esp[i1][m2])))
              {
                u[i1][m2] = interpEx(u, ell, n0, interp_points, i1, m2);
              }
            if((m2 <= n) && (isnan(esp[i][m2])))
              {
                u[i][m2] = interpEx(u, ell, n0, interp_points, i, m2);
              }
            if((i2 <= n) && (m2 <= n) && (isnan(esp[i2][m2])))
              {
                u[i2][m2] = interpEx(u, ell, n0, interp_points, i2, m2);
              }
            if((i2 <= n) && (isnan(esp[i2][m])))
              {
                u[i2][m] = interpEx(u, ell, n0, interp_points, i2, m);
              }
            if((i2 <= n) && (m1 >= 0) && (isnan(esp[i2][m1])))
              {
                u[i2][m1] = interpEx(u, ell, n0, interp_points, i2, m1);
              }
            if((m1 >= 0) && (isnan(esp[i][m1])))
              {
                u[i][m1] = interpEx(u, ell, n0, interp_points, i, m1);
              }
            if((i1 >= 0) && (m1 >= 0) && (isnan(esp[i1][m1])))
              {
                u[i1][m1] = interpEx(u, ell, n0, interp_points, i1, m1);
              }/*fim do primeiro nivel de refinamento*/
                    
                   
        
            if (ell>inf[0])
              {
                i11= i-(1<<(ell-inf[0]-1));
                i21= i+(1<<(ell-inf[0]-1));
                m11= m-(1<<(ell-inf[0]-1));
                m21= m+(1<<(ell-inf[0]-1));
                 
                if ((i11 >= 0) && (isnan(esp[i11][m])))
                  {
                    u[i11][m] = interpEx(u, ell, n0, interp_points, i11, m);
                  }
                if((i11 >= 0) && (m21 <= n) && (isnan(esp[i11][m21])))
                  {
                    u[i11][m21] = interpEx(u, ell, n0, interp_points, i11, m21);
                  }
                if((m21 <= n) && (isnan(esp[i][m21])))
                  {
                    u[i][m21] = interpEx(u, ell, n0, interp_points, i, m21);
                  }
                if((i21 <= n) && (m21 <= n) && (isnan(esp[i21][m21])))
                  {
                    u[i21][m21] = interpEx(u, ell, n0, interp_points, i21, m21);
                  }
                if((i21 <= n) && (isnan(esp[i21][m])))
                  {
                    u[i21][m] = interpEx(u, ell, n0, interp_points, i21, m);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interpEx(u, ell, n0, interp_points, i11, m11);
                  }
                if((m11 >= 0) && (isnan(esp[i][m11])))
                  {
                    u[i][m11] = interpEx(u, ell, n0, interp_points, i, m11);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interpEx(u, ell, n0, interp_points, i11, m11);
                  }
              }/*fin do segundo nivel de refinamento*/
          }/*fim dos ifs*/
                
      }/*fim do for */
         
}


/*=======================================================================*/
/* Recall that Ey is anti symmetric in the x direction 
  and symmetric in the y direction */
/*=======================================================================*/

double interpEy(double** u, int ell, int n0, int interp_points, int k, int m)
{
  int n = (n0) * (1 << ell);
  double v[interp_points];
  double h[interp_points];
  int ix[interp_points];
  int iy[interp_points];
  double r;
  int spc;
  int inf[3];
  int j;
  int level;
  int fx;
  int fy;
  find_level(inf, k, m, ell);
  level=inf[0];
  fx=inf[1];
  fy=inf[2];
  spc= 1 << (ell - level + 1);


  /*memset(v, 0, sizeof(v));
    memset(h, 0, sizeof(h));
    memset(ix, 0, sizeof(ix));
    memset(iy, 0, sizeof(iy));*/
  /*
   * The (interior) points to interpolate from are the following:
   *
   *
   * if interp_points == 4:    (need to redefine it close to the boundaries)
   *
   *    ix[0] = k - 3*spc/2;
   *    ix[1] = k - spc/2;
   *    ix[2] = k + spc/2;
   *    ix[3] = k + 3*spc/2;
   *
   * The following loop handles the general case (interp_points even and
   * greater than two).
   */

        
  if ((fx==1) && (fy==0))     /* interpolaĆ§Ć£o segundo o eixo dos XX */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }

      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[3];
          h[0] =  1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          ix[0] = ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de E\n");
              exit(-1);
            }
          if (isnan(u[ix[j]][m]))
            {
              v[j] = interpEy(u, ell, n0, interp_points, ix[j], m);
            }
          else
            {
              v[j] = u[ix[j]][m];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j]; /*faz as contas*/
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx */


  else if ((fx==0) && (fy==1))     /* interpolaĆ§Ć£o segundo o eixo dos YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          iy[interp_points/2-j] = m - (2*j-1)*spc/2;
          iy[interp_points/2-1+j] = m + (2*j-1)*spc/2;
        }

      if (iy[0] < 0) /* front boundary */
        {
          iy[0] = -iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[3];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (iy[3] > n) /* back boundary */
        {
          iy[0] = iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (iy[j] < 0 || iy[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, iy[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
          if (isnan(u[k][iy[j]]))
            {
              v[j] = interpEy(u, ell, n0, interp_points, k, iy[j]);                 
            }
          else
            {
              v[j] = u[k][iy[j]];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos yy*/

  else if ((fx==1) && (fy==1))   /* interpolaĆ§Ć£o segundo o eixo dos XX e YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }

      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[3];
          h[0] =  1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          ix[0] = ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
                
          v[j] = interpEy(u, ell, n0, interp_points, ix[j], m); 
                 
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
           
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx e dos yy */
  else
    { printf("Estou na malha mais grossa. o valor de fx e %d e de fy e %d\n", fx, fy);
      exit(-1);
    }
}
/*=======================================================================*/


void spr_data_to_sparseEy(double **u, int ell, int n0, int interp_points, double eps)
{
  int n = (n0) * (1<<ell);
  int spc = 1;
  int i,j,m;

  /* run through the data, from the finest to the coarsest detail levels */
  for (j=ell; j >= 1; j--)
    {
      spc <<= 1; /* the separation between points doubles with each level */
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpEy(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      //printf("a fiferenca e %g\n", u[i][m] - interp(u, ell, n0, interp_points, i ,m, inf));
                      //scanf("%d",&p); 
                      u[i][m] = NAN;
                           
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas linhas*/


      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpEy(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpEy(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/


    }/*termina o ciclo em j*/
}

/*=======================================================================*/

void spr_sparse_to_dataEy(double **u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1<<ell);
  int spc =  1 << (ell);
  int i,j,m;

  /* run through the data, from the finest to the coarsest detail levels */
  for (j=1; j <= ell; j++)
    {
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (isnan(u[i][m]))
                u[i][m] = interpEy(u, ell, n0, interp_points, i ,m);
              //printf("OK\n");
            }
        }/*terminou a interpolaĆ§Ć£o nas linhas*/
                  
      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interpEy(u, ell, n0, interp_points, i ,m);
                         
            }                 
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interpEy(u, ell, n0, interp_points, i ,m);
                           
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/
      spc >>= 1; /* the separation between points doubles with each level */

    }/*termina o ciclo em j*/
}

/*=======================================================================*/
void derivXEy(double** u, int ell, int n0, int interp_points, double dx, double **du,  double** v)
{

  int n = (n0) * (1<<ell);
  int i, m;

  for (i=0; i <=n; i++)
    for (m=0; m <=n; m++)   
      { if (isnan(v[i][m]))
          du[i][m] = NAN;
        else
          du[i][m]=spr_derivXEy(u, ell, n0, interp_points, dx, i, m);
      }

}
/*=======================================================================*/

double spr_derivXEy(double** u, int ell, int n0, int interp_points, double dx, int k, int m)
{
  int n = (n0) * (1 << ell);
  int ix[interp_points];        /* the positions of the samples needed to compute the derivative */
  double h[interp_points];      /* the weights needed to compute the derivative */
  double v[interp_points];      /* the function values needed to compute the derivative */
  int i;
  double du;


  
  if (k==0)
    {
      ix[0] = 2;
      ix[1] = 1;
      ix[2] = 1;
      ix[3] = 2;
      h[0] = -1./12.;
      h[1] = 2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (k==1)
    {
      ix[0] = 1;
      ix[1] = 0;
      ix[2] = 2;
      ix[3] = 3;
      h[0] = -1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (k==n-1)
    {
      ix[0] = n-3;
      ix[1] = n-2;
      ix[2] = n;
      ix[3] = n-1;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = 1./12.;
    }
  else if (k==n)
    {
      ix[0] = n-2;
      ix[1] = n-1;
      ix[2] = n-1;
      ix[3] = n-2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =-2./3.;
      h[3] = 1./12.;
    }
  else
    {
      ix[0] = k-2;
      ix[1] = k-1;
      ix[2] = k+1;
      ix[3] = k+2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] =-1./12.;
    }

  

  for (i=0; i <interp_points; i++)
    {
      if (ix[i] < 0 || ix[i] > n)
        {
          fprintf(stderr, "ERROR: spr_deriv(): ix[%d]=%d\n", i, ix[i]);
          printf("Estou na derivada de  Ex\n");
          exit(-1);
        }
      v[i] = u[ix[i]][m];
      if (isnan(v[i]))
        {
          v[i] = interpEy(u, ell, n0, interp_points, ix[i], m);         /* interpolate u[] at ix[i] */  
        }
    }
  /* compute the derivative */
  du = 0;
  for (i=0; i < interp_points; i++) du += v[i] * h[i];
  du /= (1*dx);
  return du;
}



/*=======================================================================*/

void spr_refineEy(double** u, int ell, int n0, int interp_points)
{
  //printf("estou aqui\n");
  int n = (n0) * (1 << ell);
  int i, m, i1 ,i2, i11, i21, m1, m2, m11, m21;
  int inf[3];
  double esp[n+1][n+1];
        
  //printf("estou aqui\n");

  for (i=0; i <=n; i++)
    for (m=0; m <=n; m++)
      esp[i][m]= u[i][m];

  for (m=0; m <= n;m++ )
    for (i=0; i <= n;i++ )
      {
        find_level(inf, i, m, ell);
        if ((!isnan(esp[i][m])) && (inf[0]>0))
          {
            i1= i-(1<<(ell-inf[0]));
            i2= i+(1<<(ell-inf[0]));                
            m1= m-(1<<(ell-inf[0]));
            m2= m+(1<<(ell-inf[0]));
            // printf("o valor de i1 e %d de i2 e %d de m1 e %d e de m2 e %d\n", i1, i2, m1, m2);
                 
            if ((i1 >= 0) && (isnan(esp[i1][m])))
              {
                u[i1][m] = interpEy(u, ell, n0, interp_points, i1, m);
              }
            if((i1 >= 0)  && (m2  <= n) && (isnan(esp[i1][m2])))
              {
                u[i1][m2] = interpEy(u, ell, n0, interp_points, i1, m2);
              }
            if((m2 <= n) && (isnan(esp[i][m2])))
              {
                u[i][m2] = interpEy(u, ell, n0, interp_points, i, m2);
              }
            if((i2 <= n) && (m2 <= n) && (isnan(esp[i2][m2])))
              {
                u[i2][m2] = interpEy(u, ell, n0, interp_points, i2, m2);
              }
            if((i2 <= n) && (isnan(esp[i2][m])))
              {
                u[i2][m] = interpEy(u, ell, n0, interp_points, i2, m);
              }
            if((i2 <= n) && (m1 >= 0) && (isnan(esp[i2][m1])))
              {
                u[i2][m1] = interpEy(u, ell, n0, interp_points, i2, m1);
              }
            if((m1 >= 0) && (isnan(esp[i][m1])))
              {
                u[i][m1] = interpEy(u, ell, n0, interp_points, i, m1);
              }
            if((i1 >= 0) && (m1 >= 0) && (isnan(esp[i1][m1])))
              {
                u[i1][m1] = interpEy(u, ell, n0, interp_points, i1, m1);
              }/*fim do primeiro nivel de refinamento*/
                    
                   
        
            if (ell>inf[0])
              {
                i11= i-(1<<(ell-inf[0]-1));
                i21= i+(1<<(ell-inf[0]-1));
                m11= m-(1<<(ell-inf[0]-1));
                m21= m+(1<<(ell-inf[0]-1));
                 
                if ((i11 >= 0) && (isnan(esp[i11][m])))
                  {
                    u[i11][m] = interpEy(u, ell, n0, interp_points, i11, m);
                  }
                if((i11 >= 0) && (m21 <= n) && (isnan(esp[i11][m21])))
                  {
                    u[i11][m21] = interpEy(u, ell, n0, interp_points, i11, m21);
                  }
                if((m21 <= n) && (isnan(esp[i][m21])))
                  {
                    u[i][m21] = interpEy(u, ell, n0, interp_points, i, m21);
                  }
                if((i21 <= n) && (m21 <= n) && (isnan(esp[i21][m21])))
                  {
                    u[i21][m21] = interpEy(u, ell, n0, interp_points, i21, m21);
                  }
                if((i21 <= n) && (isnan(esp[i21][m])))
                  {
                    u[i21][m] = interpEy(u, ell, n0, interp_points, i21, m);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interpEy(u, ell, n0, interp_points, i11, m11);
                  }
                if((m11 >= 0) && (isnan(esp[i][m11])))
                  {
                    u[i][m11] = interpEy(u, ell, n0, interp_points, i, m11);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interpEy(u, ell, n0, interp_points, i11, m11);
                  }
              }/*fin do segundo nivel de refinamento*/
          }/*fim dos ifs*/
                
      }/*fim do for */
         
}


/*=======================================================================*/

/* LOCAL function
 *
 * This function does the basic interpolation of Hz: it computes the value u[k],
 * given the sparse representation stored in u[0..n]. It is tacitly assumed
 * that u[k] is not in the sparse decomposition, that is, u[k] is NAN. The
 * number of points used to interpolate is interp_points, and the spacing at
 * the next coarsest level is spc.
 */
double interpH(double **u, int ell, int n0, int interp_points, int k, int m)
{
  int n = (n0) * (1 << ell);
  double v[interp_points];
  double h[interp_points];
  int ix[interp_points];
  int iy[interp_points];
  double r;
  int j;
  int inf[3];
  int level, spc;
  int fx;
  int fy;
  find_level(inf, k, m, ell);
  fx=inf[1];
  fy=inf[2];
  level=inf[0];
  spc = 1 << (ell - level + 1);

  /*memset(v, 0, sizeof(v));
    memset(h, 0, sizeof(h));
    memset(ix, 0, sizeof(ix));
    memset(iy, 0, sizeof(iy));*/
  /*
   * The (interior) points to interpolate from are the following:
   *
   *
   * if interp_points == 4:    (need to redefine it close to the boundaries)
   *
   *    ix[0] = k - 3*spc/2;
   *    ix[1] = k - spc/2;
   *    ix[2] = k + spc/2;
   *    ix[3] = k + 3*spc/2;
   *
   * The following loop handles the general case (interp_points even and
   * greater than two).
   */

        
  if ((fx==1) && (fy==0))     /* interpolaĆ§Ć£o segundo o eixo dos XX */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }
      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[3];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          ix[0] = ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de H\n");
              exit(-1);
            }
          if (isnan(u[ix[j]][m]))
            {
              v[j] = interpH(u, ell, n0, interp_points, ix[j], m);
            }
          else
            {
              v[j] = u[ix[j]][m];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j]; /*faz as contas*/
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx */


  else if ((fx==0) && (fy==1))     /* interpolaĆ§Ć£o segundo o eixo dos YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          iy[interp_points/2-j] = m - (2*j-1)*spc/2;
          iy[interp_points/2-1+j] = m + (2*j-1)*spc/2;
        }

      if (iy[0] < 0) /* front boundary */
        {
          iy[0] = -iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[3];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (iy[3] > n) /* back boundary */
        {
          iy[0] = iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (iy[j] < 0 || iy[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, iy[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
          if (isnan(u[k][iy[j]]))
            {
              v[j] = interpH(u, ell, n0, interp_points, k, iy[j]); 
            }
          else
            {
              v[j] = u[k][iy[j]];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos yy*/

  else  if ((fx==1) && (fy==1)) /*interpolaĆ§Ć£o segundo o eixo dos XX e YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }

      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[3];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          ix[0] = ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
                
          v[j] = interp(u, ell, n0, interp_points,  ix[j], m); 
                
                
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx e dos yy */
  else 
    { printf("Estou na malha mais grossa. o valor de fx e %d e de fy e %d e o nivel e %d\n", fx, fy,level);
      exit(-1);
    }

}
/*=======================================================================*/

void spr_data_to_sparseH(double** u, int ell, int n0, int interp_points, double eps)
{
  int n = (n0) * (1<<ell);
  int spc = 1;
  int i,j,m;
        

  /* run through the data, from the finest to the coarsest detail levels */
  for (j=ell; j >= 1; j--)
    {
      spc <<= 1; /* the separation between points doubles with each level */
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpH(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas linhas*/


      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpH(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interpH(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/


    }/*termina o ciclo em j*/
}

/*=======================================================================*/

void spr_sparse_to_dataH(double **u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1<<ell);
  int spc =  1 << (ell);
  int i,j,m;
        

  /* run through the data, from the finest to the coarsest detail levels */
  for (j=1; j <= ell; j++)
    {
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (isnan(u[i][m]))
                u[i][m] = interpH(u, ell, n0, interp_points, i ,m);
            }
        }/*terminou a interpolaĆ§Ć£o nas linhas*/
                  
      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interpH(u, ell, n0, interp_points, i ,m);
                         
            }                 
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interpH(u, ell, n0, interp_points, i ,m);
                           
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/
      spc >>= 1; /* the separation between points doubles with each level */

    }/*termina o ciclo em j*/
}

/*=======================================================================*/

double spr_derivXH(double** u, int ell, int n0, int interp_points, double dx, int k, int m)
{
  int n = (n0) * (1 << ell);
  int ix[interp_points];        /* the positions of the samples needed to compute the derivative */
  double h[interp_points];      /* the weights needed to compute the derivative */
  double v[interp_points];      /* the function values needed to compute the derivative */
  int i;
  double du;
  

  if (k==0)
    {
      ix[0] = 2;
      ix[1] = 1;
      ix[2] = 1;
      ix[3] = 2;
      h[0] = 1./12.;
      h[1] = -2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (k==1)
    {
      ix[0] = 1;
      ix[1] = 0;
      ix[2] = 2;
      ix[3] = 3;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (k==n-1)
    {
      ix[0] = n-3;
      ix[1] = n-2;
      ix[2] = n;
      ix[3] = n-1;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (k==n)
    {
      ix[0] = n-2;
      ix[1] = n-1;
      ix[2] = n-1;
      ix[3] = n-2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =2./3.;
      h[3] =-1./12.;
    }
  else
    {
      ix[0] = k-2;
      ix[1] = k-1;
      ix[2] = k+1;
      ix[3] = k+2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] =-1./12.;
    }

  

  for (i=0; i <interp_points; i++)
    {
      if (ix[i] < 0 || ix[i] > n)
        {
          fprintf(stderr, "ERROR: spr_deriv(): ix[%d]=%d\n", i, ix[i]);
          printf("Estou na derivada de  Ex\n");
          exit(-1);
        }
      v[i] = u[ix[i]][m];
      if (isnan(v[i]))
        {
          v[i] = interpH(u, ell, n0, interp_points, ix[i], m);  /* interpolate u[] at ix[i] */  
        }
    }
  /* compute the derivative */
  du = 0;
  for (i=0; i < interp_points; i++) du += v[i] * h[i];
  du /= (1*dx);
  return du;
}



/*=======================================================================*/

double spr_derivYH(double** u, int ell, int n0, int interp_points, double dx, int k, int m)
{
  int n = (n0) * (1 << ell);
  int ix[interp_points];        /* the positions of the samples needed to compute the derivative */
  double h[interp_points];      /* the weights needed to compute the derivative */
  double v[interp_points];      /* the function values needed to compute the derivative */
  int i;
  double du;

  if (m==0)
    {
      ix[0] = 2;
      ix[1] = 1;
      ix[2] = 1;
      ix[3] = 2;
      h[0] = 1./12.;
      h[1] = -2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (m==1)
    {
      ix[0] = 1;
      ix[1] = 0;
      ix[2] = 2;
      ix[3] = 3;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (m==n-1)
    {
      ix[0] = n-3;
      ix[1] = n-2;
      ix[2] = n;
      ix[3] = n-1;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (m==n)
    {
      ix[0] = n-2;
      ix[1] = n-1;
      ix[2] = n-1;
      ix[3] = n-2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =2./3.;
      h[3] =-1./12.;
    }
  else
    {
      ix[0] = m-2;
      ix[1] = m-1;
      ix[2] = m+1;
      ix[3] = m+2;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] =-1./12.;
    }

  

  for (i=0; i <interp_points; i++)
    {
      if (ix[i] < 0 || ix[i] > n)
        {
          fprintf(stderr, "ERROR: spr_deriv(): ix[%d]=%d\n", i, ix[i]);
          printf("Estou na derivada de  Ex\n");
          exit(-1);
        }
      v[i] = u[k][ix[i]];
      if (isnan(v[i]))
        {
          v[i] = interpH(u, ell, n0, interp_points, k, ix[i]);  /* interpolate u[] at ix[i] */  
        }
    }
  /* compute the derivative */
  du = 0;
  for (i=0; i < interp_points; i++) du += v[i] * h[i];
  du /= (1*dx);
  return du;
}

/*=======================================================================*/

void derivXH(
             double** u, /* the sparse representation u[0..n-1] */
             int ell, int n0, int interp_points,
             double dx, /* separation between data points */
             double **du,   /* the sparse representation of derivative */
             double** v)
{

  int n = (n0) * (1<<ell);
  int i, m;

  for (i=1; i <=n-1; i++)
    for (m=0; m <=n; m++)   
      { if (isnan(v[i][m]))
          du[i][m] = NAN;
        else
          du[i][m]=spr_derivXH(u, ell, n0, interp_points, dx, i, m);
      }
  for (m=0; m <=n; m++)
    du[0][m]=du[n][m]=0.0;

}
/*=======================================================================*/
void derivYH(
             double** u, /* the sparse representation u[0..n-1] */
             int ell, int n0, int interp_points,
             double dx, /* separation between data points */
             double **du,   /* the sparse representation of derivative */
             double** v)
{

  int n = (n0) * (1<<ell);
  int i, m;

  for (i=0; i <=n; i++)
    for (m=1; m <=n-1; m++) 
      { if (isnan(v[i][m]))
          du[i][m] = NAN;
        else
          du[i][m]=spr_derivYH(u, ell, n0, interp_points, dx, i, m);
      }
  for (i=0; i <=n; i++)
    du[i][0]=du[i][n]=0.0;

}
/*=======================================================================*/

void spr_refine(double** u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1 << ell);
  int i, m, i1 ,i2, i11, i21, m1, m2, m11, m21;
  int inf[3];
  double esp[n+1][n+1];
  for (i=0; i <=n; i++)
    for (m=0; m <=n; m++)
      esp[i][m]= u[i][m];
                

  for (m=1; m <= n-1;m++ )
    for (i=1; i <= n-1;i++ )
      {
        if (!isnan(esp[i][m]))
          {
            find_level(inf, i, m, ell);
            i1= i-(1<<(ell-inf[0]));
            i2= i+(1<<(ell-inf[0]));                
            m1= m-(1<<(ell-inf[0]));
            m2= m+(1<<(ell-inf[0]));
                 
            if ((i1 >= 0) && (isnan(esp[i1][m])))
              {
                u[i1][m] = interp(u, ell, n0, interp_points, i1, m);
              }
            if((i1 >= 0) && (m2 <= n) && (isnan(esp[i1][m2])))
              {
                u[i1][m2] = interp(u, ell, n0, interp_points, i1, m2);
              }
            if((m2 >= 0) && (isnan(esp[i][m2])))
              {
                u[i][m2] = interp(u, ell, n0, interp_points, i, m2);
              }
            if((i2 <= n) && (m2 <= n) && (isnan(esp[i2][m2])))
              {
                u[i2][m2] = interp(u, ell, n0, interp_points, i2, m2);
              }
            if((i2 <= n) && (isnan(esp[i2][m])))
              {
                u[i2][m] = interp(u, ell, n0, interp_points, i2, m);
              }
            if((i2 <= n) && (m1 >= 0) && (isnan(esp[i2][m1])))
              {
                u[i2][m1] = interp(u, ell, n0, interp_points, i2, m1);
              }
            if((m1 >= 0) && (isnan(esp[i][m1])))
              {
                u[i][m1] = interp(u, ell, n0, interp_points, i, m1);
              }
            if((i1 >= 0) && (m1 >= 0) && (isnan(esp[i1][m1])))
              {
                u[i1][m1] = interp(u, ell, n0, interp_points, i1, m1);
              }/*fim do primeiro nivel de refinamento*/
                    
                   
        
            if (ell-inf[0])
              {
                i11= i-(1<<(ell-inf[0]-1));
                i21= i+(1<<(ell-inf[0]-1));
                m11= m-(1<<(ell-inf[0]-1));
                m21= m+(1<<(ell-inf[0]-1));

 
                 
                if ((i11 >= 0) && (isnan(esp[i11][m])))
                  {
                    u[i11][m] = interp(u, ell, n0, interp_points, i11, m);
                  }
                if((i11 >= 0) && (m21 <= n) && (isnan(esp[i11][m21])))
                  {
                    u[i11][m21] = interp(u, ell, n0, interp_points, i11, m21);
                  }
                if((m21 <= n) && (isnan(esp[i][m21])))
                  {
                    u[i][m21] = interp(u, ell, n0, interp_points, i, m21);
                  }
                if((i21 <= n) && (m21 <= n) && (isnan(esp[i21][m21])))
                  {
                    u[i21][m21] = interp(u, ell, n0, interp_points, i21, m21);
                  }
                if((i21 <= n) && (isnan(esp[i21][m])))
                  {
                    u[i21][m] = interp(u, ell, n0, interp_points, i21, m);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interp(u, ell, n0, interp_points, i11, m11);
                  }
                if((m11 >= 0) && (isnan(esp[i][m11])))
                  {
                    u[i][m11] = interp(u, ell, n0, interp_points, i, m11);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interp(u, ell, n0, interp_points, i11, m11);
                  }
              }/*fin do segundo nivel de refinamento*/
          }/*fim dos ifs*/
                
      }/*fim do for */
         
}

/*=======================================================================*/

//extern void** dvector(int,int);

void spr_refineH(double** u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1 << ell);
  int i, m, i1 ,i2, i11, i21, m1, m2, m11, m21;
  int inf[3];
  double esp[n+1][n+1];
  for (i=0; i <=n; i++)
    for (m=0; m <=n; m++)
      esp[i][m]= u[i][m];
                

  for (m=0; m <= n;m++ )
    for (i=0; i <= n;i++ )
      {
        find_level(inf, i, m, ell);
        if ((!isnan(esp[i][m])) && (inf[0]>0))
          {
            i1= i-(1<<(ell-inf[0]));
            i2= i+(1<<(ell-inf[0]));                
            m1= m-(1<<(ell-inf[0]));
            m2= m+(1<<(ell-inf[0]));
                 
            if ((i1 >= 0) && (isnan(esp[i1][m])))
              {
                u[i1][m] = interpH(u, ell, n0, interp_points, i1, m);
              }
            if((i1 >= 0)  && (m2  <= n) && (isnan(esp[i1][m2])))
              {
                u[i1][m2] = interpH(u, ell, n0, interp_points, i1, m2);
              }
            if((m2 <= n) && (isnan(esp[i][m2])))
              {
                u[i][m2] = interpH(u, ell, n0, interp_points, i, m2);
              }
            if((i2 <= n) && (m2 <= n) && (isnan(esp[i2][m2])))
              {
                u[i2][m2] = interpH(u, ell, n0, interp_points, i2, m2);
              }
            if((i2 <= n) && (isnan(esp[i2][m])))
              {
                u[i2][m] = interpH(u, ell, n0, interp_points, i2, m);
              }
            if((i2 <= n) && (m1 >= 0) && (isnan(esp[i2][m1])))
              {
                u[i2][m1] = interpH(u, ell, n0, interp_points, i2, m1);
              }
            if((m1 >= 0) && (isnan(esp[i][m1])))
              {
                u[i][m1] = interpH(u, ell, n0, interp_points, i, m1);
              }
            if((i1 >= 0) && (m1 >= 0) && (isnan(esp[i1][m1])))
              {
                u[i1][m1] = interpH(u, ell, n0, interp_points, i1, m1);
              }/*fim do primeiro nivel de refinamento*/
                    
                   
        
            if (ell>inf[0])
              {
                i11= i-(1<<(ell-inf[0]-1));
                i21= i+(1<<(ell-inf[0]-1));
                m11= m-(1<<(ell-inf[0]-1));
                m21= m+(1<<(ell-inf[0]-1));
                 
                if ((i11 >= 0) && (isnan(esp[i11][m])))
                  {
                    u[i11][m] = interpH(u, ell, n0, interp_points, i11, m);
                  }
                if((i11 >= 0) && (m21 <= n) && (isnan(esp[i11][m21])))
                  {
                    u[i11][m21] = interpH(u, ell, n0, interp_points, i11, m21);
                  }
                if((m21 <= n) && (isnan(esp[i][m21])))
                  {
                    u[i][m21] = interpH(u, ell, n0, interp_points, i, m21);
                  }
                if((i21 <= n) && (m21 <= n) && (isnan(esp[i21][m21])))
                  {
                    u[i21][m21] = interpH(u, ell, n0, interp_points, i21, m21);
                  }
                if((i21 <= n) && (isnan(esp[i21][m])))
                  {
                    u[i21][m] = interpH(u, ell, n0, interp_points, i21, m);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interpH(u, ell, n0, interp_points, i11, m11);
                  }
                if((m11 >= 0) && (isnan(esp[i][m11])))
                  {
                    u[i][m11] = interpH(u, ell, n0, interp_points, i, m11);
                  }
                if((i11 >= 0) && (m11 >= 0) && (isnan(esp[i11][m11])))
                  {
                    u[i11][m11] = interpH(u, ell, n0, interp_points, i11, m11);
                  }
              }/*fin do segundo nivel de refinamento*/
          }/*fim dos ifs*/
                
      }/*fim do for */
         
}



/*=======================================================================*/
int find_level_max(double** u, int ell, int n0, int ell_atual)
{
  int n = (n0) * (1 << ell);
  int i, m, j;
  int spc = 1 << (ell - ell_atual);
  for (j=ell_atual; j >= 1; j--)
    {
      spc <<= 1; /* the separation between points doubles with each level */
      for (m=0; m <=n; m+=spc)
        for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
          if (!isnan(u[i][m])) return j;      
      /*terminou de percorrer as linhas*/

    
      for (i=0; i <=n; i+=spc)
        for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
          if (!isnan(u[i][m])) return j;
      /*terminou de percorrer as colunas*/ 

      for (i=spc/2; i <n; i+=spc)
        for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
          if (!isnan(u[i][m]))    return j;
      /*terminou de percorrer os pontos do meio*/        
                    
    }/*termina o ciclo em j*/
           
  return 0;
}

/*=======================================================================================*/

/*=======================================================================*/

/* LOCAL function
 *
 * This function does the basic interpolation of a function with antisymmetric boundary conditions
 * in all sides: it computes the value u[k],
 * given the sparse representation stored in u[0..n]. It is tacitly assumed
 * that u[k] is not in the sparse decomposition, that is, u[k] is NAN. The
 * number of points used to interpolate is interp_points, and the spacing at
 * the next coarsest level is spc.
 */
double interp(double** u, int ell, int n0, int interp_points, int k, int m)
{
  int n = (n0) * (1 << ell);
  double v[interp_points];
  double h[interp_points];
  int ix[interp_points];
  int iy[interp_points];
  double r;
  int spc;
  int inf[3];
  int j;
  int level;
  int fx;
  int fy;
  find_level(inf, k, m, ell);
  level=inf[0];
  fx=inf[1];
  fy=inf[2];
  spc= 1 << (ell - level + 1);


  /*memset(v, 0, sizeof(v));
    memset(h, 0, sizeof(h));
    memset(ix, 0, sizeof(ix));
    memset(iy, 0, sizeof(iy));*/
  /*
   * The (interior) points to interpolate from are the following:
   *
   *
   * if interp_points == 4:    (need to redefine it close to the boundaries)
   *
   *    ix[0] = k - 3*spc/2;
   *    ix[1] = k - spc/2;
   *    ix[2] = k + spc/2;
   *    ix[3] = k + 3*spc/2;
   *
   * The following loop handles the general case (interp_points even and
   * greater than two).
   */

        
  if ((fx==1) && (fy==0))     /* interpolaĆ§Ć£o segundo o eixo dos XX */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }

      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[3];
          h[0] =  1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          ix[0] = ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de E\n");
              exit(-1);
            }
          if (isnan(u[ix[j]][m]))
            {
              v[j] = interp(u, ell, n0, interp_points, ix[j], m);
            }
          else
            {
              v[j] = u[ix[j]][m];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j]; /*faz as contas*/
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx */


  else if ((fx==0) && (fy==1))     /* interpolaĆ§Ć£o segundo o eixo dos YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          iy[interp_points/2-j] = m - (2*j-1)*spc/2;
          iy[interp_points/2-1+j] = m + (2*j-1)*spc/2;
        }

      if (iy[0] < 0) /* left boundary */
        {
          iy[0] = -iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[3];
          h[0] =  1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (iy[3] > n) /* right boundary */
        {
          iy[0] = iy[0];
          iy[1] = iy[1];
          iy[2] = iy[2];
          iy[3] = iy[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (iy[j] < 0 || iy[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, iy[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
          if (isnan(u[k][iy[j]]))
            {
              v[j] = interp(u, ell, n0, interp_points, k, iy[j]);           
            }
          else
            {
              v[j] = u[k][iy[j]];
            }
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos yy*/

  else if ((fx==1) && (fy==1))   /* interpolaĆ§Ć£o segundo o eixo dos XX e YY */
    {

      for (j=1; j <= interp_points/2; j++)
        {
          ix[interp_points/2-j] = k - (2*j-1)*spc/2;
          ix[interp_points/2-1+j] = k + (2*j-1)*spc/2;
        }

      if (ix[0] < 0) /* left boundary */
        {
          ix[0] = -ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[3];
          h[0] =  1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  -1/16.0;
        }
      else if (ix[3] > n) /* right boundary */
        {
          ix[0] = ix[0];
          ix[1] = ix[1];
          ix[2] = ix[2];
          ix[3] = ix[1];
          h[0] =  -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] =  1/16.0;
        }
      else /* interior */
        {
          h[0] = -1/16.0;
          h[1] =  9/16.0;
          h[2] =  9/16.0;
          h[3] = -1/16.0;
        }

      for (j=0; j < interp_points; j++)
        {
          if (ix[j] < 0 || ix[j] > n)
            {
              fprintf(stderr, "ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              printf("Estou na interpolacao para a malha de  E\n");
              exit(-1);
            }
                
          v[j] = interp(u, ell, n0, interp_points, ix[j], m); 
                 
        }

      for (r=j=0; j < interp_points; j++) r += h[j] * v[j];
      return r;
           
    } /*fim da interpolaĆ§Ć£o segundo o eixo dos xx e dos yy */
  else
    { printf("Estou na malha mais grossa. o valor de fx e %d e de fy e %d\n", fx, fy);
      exit(-1);
    }
}
/*=======================================================================*/

/*
 * PUBLIC function
 *
 * Converts a data vector of n elements to a sparse representation.
 * It sets a data value u[i] in u[0..n-1] to NAN if the absolute value of the
 * corresponding wavelet coefficient d[i] is smaller than epsilon.
 */
void spr_data_to_sparse(double **u, int ell, int n0, int interp_points, double eps)
{
  int n = (n0) * (1<<ell);
  int spc = 1;
  int i,j,m;

  /* run through the data, from the finest to the coarsest detail levels */
  for (j=ell; j >= 1; j--)
    {
      spc <<= 1; /* the separation between points doubles with each level */
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interp(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      //printf("a fiferenca e %g\n", u[i][m] - interp(u, ell, n0, interp_points, i ,m, inf));
                      //scanf("%d",&p); 
                      u[i][m] = NAN;
                           
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas linhas*/


      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interp(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (!isnan(u[i][m]))
                {
                  if (fabs(u[i][m] - interp(u, ell, n0, interp_points, i ,m)) < eps)
                    {
                      u[i][m] = NAN;
                    }
                }
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/


    }/*termina o ciclo em j*/
}

/*=======================================================================*/

void spr_sparse_to_data(double **u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1<<ell);
  int spc =  1 << (ell);
  int i,j,m;

  /* run through the data, from the finest to the coarsest detail levels */
  for (j=1; j <= ell; j++)
    {
      for (m=0; m <=n; m+=spc)
        { 
          for (i=spc/2; i < n; i+=spc) /* run through all points at level j e da linha m */
            {
              if (isnan(u[i][m]))
                u[i][m] = interp(u, ell, n0, interp_points, i ,m);
              //printf("OK\n");
            }
        }/*terminou a interpolaĆ§Ć£o nas linhas*/
                  
      for (i=0; i <=n; i+=spc)
        {
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interp(u, ell, n0, interp_points, i ,m);
                         
            }                 
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      for (i=spc/2; i <n; i+=spc)
        { 
          for (m=spc/2; m < n; m+=spc) /* run through all points at level j e da linha m coluna i*/
            {
              if (isnan(u[i][m]))
                u[i][m] = interp(u, ell, n0, interp_points, i ,m);
                           
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/
      spc >>= 1; /* the separation between points doubles with each level */

    }/*termina o ciclo em j*/
}