/* Last edited on 2008-05-21 18:13:27 by anamaria */ 

/*
 * 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 "Int2Dlocal.h"
int cont;

 

/*
 * 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).
   */

   cont++;
        
  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;

 
  int escala; 

 escala=local_scaley(u, k, m,  ell, n0);
 
	           

/* if (k==3)
 {    fprintf(stderr, "escala %i\n", escala);
 }*/
  
  if (m==0)
    {
      ix[0] = 2*escala;
      ix[1] = escala;
      ix[2] = escala;
      ix[3] = 2*escala;
      h[0] = -1./12.;
      h[1] = +2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if ((m<=(2*escala))&&((m+2*escala)<=n))
    {
      ix[0] = 2*escala-m;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = m+2*escala;
      h[0] = -1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  
  else if (m==n)
    {
      ix[0] = n-2*escala;
      ix[1] = n-escala;
      ix[2] = n-escala;
      ix[3] = n-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =-2./3.;
      h[3] =+1./12.;
    }
    else if ((m>(n-2*escala))&&((m-2*escala)>=0))
    {
      ix[0] = m-2*escala;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = 2*n-m-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = +1./12.;
    }
    else if ((m>(n-2*escala))&&((m-2*escala)<0))
    {
    ix[0] = 2*escala-m;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = 2*n-m-2*escala;
      h[0] = -1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = 1./12.;
    }
  else
    {
      ix[0] = m-2*escala;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = m+2*escala;
      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 /= (escala*dx);
  return du;
}

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

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=dvector(n+1, n+1);
        
                

  for (i=0; i <=n; i++)
    for (m=0; m <=n; m++)
      esp[i][m]= u[i][m];
  //printf("Estou aqui1");
  for (m=0; m <= n;m++ )
    for (i=0; i <= n;i++ )
      { //printf("Estou aqui2");
        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((i21 >= 0) && (i21 <= n) && (m11 >= 0) && (isnan(esp[i21][m11])))
                  {
                    u[i21][m11] = interpEx(u, ell, n0, interp_points, i21, m11);
                  }
              }/*fin do segundo nivel de refinamento*/
          }/*fim dos ifs*/
                
      }/*fim do for */
    //     for (i=0; i <=n; i++)
      //    	for (m=0; m <=n; m++)
      			libera(n+1,n+1,esp);
}


/*=======================================================================*/
/* 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);

cont++;

  /*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;


  int escala; 

 escala=local_scalex(u, k, m,  ell, n0);
 
	           

/* if (k==3)
 {    fprintf(stderr, "escala %i\n", escala);
 }*/
  
  if (k==0)
    {
      ix[0] = 2*escala;
      ix[1] = escala;
      ix[2] = escala;
      ix[3] = 2*escala;
      h[0] = -1./12.;
      h[1] = +2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if ((k<=(2*escala))&&((k+2*escala)<=n))
    {
      ix[0] = 2*escala-k;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = k+2*escala;
      h[0] = -1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  
  else if (k==n)
    {
      ix[0] = n-2*escala;
      ix[1] = n-escala;
      ix[2] = n-escala;
      ix[3] = n-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =-2./3.;
      h[3] =+1./12.;
    }
    else if ((k>(n-2*escala))&&((k-2*escala)>=0))
    {
      ix[0] = k-2*escala;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = 2*n-k-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = +1./12.;
    }
    else if ((k>(n-2*escala))&&((k-2*escala)<0))
    {
    ix[0] = 2*escala-k;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = 2*n-k-2*escala;
      h[0] = -1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = 1./12.;
    }
  else
    {
      ix[0] = k-2*escala;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = k+2*escala;
      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 /= (escala*dx);
  return du;
}



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

void spr_refineEy(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=dvector(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] = 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((i21 >= 0) && (i21 <= n) && (m11 >= 0) && (isnan(esp[i21][m11])))
                  {
                    u[i21][m11] = interpEy(u, ell, n0, interp_points, i21, m11);
                  }
              }/*fin do segundo nivel de refinamento*/
          }/*fim dos ifs*/
                
      }/*fim do for */
     //     for (i=0; i <=n; i++)
       //   	for (m=0; m <=n; m++)
      			libera(n+1,n+1,esp);
}


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

/* 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;
  cont++;
  find_level(inf, k, m, ell);
  fx=inf[1];
  fy=inf[2];
  level=inf[0];
  spc = 1 << (ell - level + 1);

cont++;
  
  /*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)
            {
              printf("ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              //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)
            {
              printf("ERROR: interp(): ix[%d]=%d\n", j, iy[j]);
              //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)
            {
              printf("ERROR: interp(): ix[%d]=%d\n", j, ix[j]);
              //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;
  int escala; 

 escala=local_scalex(u, k, m,  ell, n0);
 
	           

/* if (k==3)
 {    fprintf(stderr, "escala %i\n", escala);
 }*/
  
  if (k==0)
    {
      ix[0] = 2*escala;
      ix[1] = escala;
      ix[2] = escala;
      ix[3] = 2*escala;
      h[0] = 1./12.;
      h[1] = -2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if ((k<=(2*escala))&&((k+2*escala)<=n))
    {
      ix[0] = 2*escala-k;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = k+2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  
  else if (k==n)
    {
      ix[0] = n-2*escala;
      ix[1] = n-escala;
      ix[2] = n-escala;
      ix[3] = n-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =2./3.;
      h[3] =-1./12;
    }
    else if ((k>(n-2*escala))&&((k-2*escala)>=0))
    {
      ix[0] = k-2*escala;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = 2*n-k-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
    else if ((k>(n-2*escala))&&((k-2*escala)<0))
    {
    ix[0] = 2*escala-k;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = 2*n-k-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else
    {
      ix[0] = k-2*escala;
      ix[1] = k-escala;
      ix[2] = k+escala;
      ix[3] = k+2*escala;
      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)
        {
          printf("ERROR: spr_deriv(): %i ix[%i]=%i %i\n",k, i,  ix[i], escala);
          //fprintf(stderr, "ERROR: spr_deriv(): %i ix[%i]=%i %i\n",k, i,  ix[i], escala);
          printf("Estou na derivada x de  H\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 /= (escala*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;
  int escala; 
  

 escala=local_scaley(u, k, m,  ell, n0);
 
  
  if (m==0)

    {
      ix[0] = 2*escala;
      ix[1] = escala;
      ix[2] = escala;
      ix[3] = 2*escala;
      h[0] = 1./12.;
      h[1] = -2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if ((m<=(2*escala))&&((m+2*escala)<=n))
    {
      ix[0] = 2*escala-m;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = m+2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else if (m==n)
    {
      ix[0] = n-2*escala;
      ix[1] = n-escala;
      ix[2] = n-escala;
      ix[3] = n-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] =2./3.;
      h[3] =-1./12.;
    }
        else if ((m>(n-2*escala))&&((m-2*escala)>=0))
    {
      ix[0] = m-2*escala;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = 2*n-m-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
    else if ((m>(n-2*escala))&&((m-2*escala)<0))
    {
    ix[0] = 2*escala-m;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = 2*n-m-2*escala;
      h[0] = 1./12.;
      h[1] =-2./3.;
      h[2] = 2./3.;
      h[3] = -1./12.;
    }
  else
    {
      ix[0] = m-2*escala;
      ix[1] = m-escala;
      ix[2] = m+escala;
      ix[3] = m+2*escala;
      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)
        {
          printf("ERROR: spr_deriv(): %i ix[%i]=%i %i\n",m, i,  ix[i], escala);
          //fprintf(stderr, "ERROR: spr_deriv(): %i ix[%i]=%i %i\n",m, i,  ix[i], escala);
          printf("Estou na derivada y de  H\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 /= (escala*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((i21 >= 0) && (i21 <= n) && (m11 >= 0) && (isnan(esp[i21][m11])))
                  {
                    u[i21][m11] = interp(u, ell, n0, interp_points, i21, m11);
                  }
              }/*fin do segundo nivel de refinamento*/
          }/*fim dos ifs*/
                
      }/*fim do for */
         
}

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

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=dvector(n+1, n+1);
 // printf("Estou aqui 1\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]));
                 
            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*/
                    
           //printf("estou aqui\n");        
        
               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((i21 >= 0) && (i21 <= n) && (m11 >= 0) && (isnan(esp[i21][m11])))
                  {
                    u[i21][m11] = interpH(u, ell, n0, interp_points, i21, m11);
                  }
                  } /*fin do segundo nivel de refinamento*/
                 
		//printf("estou aqui lixado\n");
          }/*fim dos ifs*/
                
      }/*fim do for */ 
  //for (i=0; i <=n; i++)
    //      	for (m=0; m <=n; m++)
      			libera(n+1,n+1,esp);
         
}



/*=======================================================================*/
void spr_coneH(double** u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1<<ell);
  int spc = 1;
  int i, i1, i2, j, m, m1, m2;
        

  /* 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]))
                {
                  
		i1= i-spc/2;
                i2= i+spc/2; 
		if ((i1 >= 0) && (isnan(u[i1][m])))
              {
                u[i1][m] = interpH(u, ell, n0, interp_points, i1, m);
              }  
	       if((i2 <= n)  && (isnan(u[i2][m])))
                  {
                    u[i2][m] = interpH(u, ell, n0, interp_points, i2, 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]))
                {
		m1= m-spc/2;
                m2= m+spc/2; 
		if ((m1 >= 0) && (isnan(u[i][m1])))
              {
                u[i][m1] = interpH(u, ell, n0, interp_points, i, m1);
              }  
	       if((m2 <= n)  && (isnan(u[i][m2])))
                  {
                    u[i][m2] = interpH(u, ell, n0, interp_points, i, m2);
                  }       
                  
                }
            }
        } /*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]))
                {
		i1= i-spc/2;
                i2= i+spc/2; 
		m1= m-spc/2;
                m2= m+spc/2; 
		   
                if((i1 >= 0) && (m2 <= n) && (isnan(u[i1][m2])))
                  {
                    u[i1][m2] = interpH(u, ell, n0, interp_points, i1, m2);
                  }
                
                if((i2 <= n) && (m2 <= n) && (isnan(u[i2][m2])))
                  {
                    u[i2][m2] = interpH(u, ell, n0, interp_points, i2, m2);
                  }
		 if((i1 >= 0) && (m2 <= n) && (isnan(u[i1][m2])))
                  {
                    u[i1][m2] = interpH(u, ell, n0, interp_points, i1, m2);
                  }
                
                if((i1 >= 0) && (m1 >= 0) && (isnan(u[i1][m1])))
                  {
                    u[i1][m1] = interpH(u, ell, n0, interp_points, i1, m1);
                  } 
                  
                }
            }
        } /*terminou a interpolaĆ§Ć£o no meio*/


    }/*termina o ciclo em j*/
}


/*=======================================================================*/
void spr_coneHSimp(double** u, int ell, int n0, int interp_points)
{
  int n = (n0) * (1<<ell);
  int spc = 1;
  int i, i1, i2, j, m, m1, m2;
        

  /* 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]))
                {
                  
		i1= i-spc/2;
                i2= i+spc/2; 
		if ((i1 >= 0) && (isnan(u[i1][m])))
              {
                u[i1][m] = interpH(u, ell, n0, interp_points, i1, m);
              }  
	       if((i2 <= n)  && (isnan(u[i2][m])))
                  {
                    u[i2][m] = interpH(u, ell, n0, interp_points, i2, 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]))
                {
		m1= m-spc/2;
                m2= m+spc/2; 
		if ((m1 >= 0) && (isnan(u[i][m1])))
              {
                u[i][m1] = interpH(u, ell, n0, interp_points, i, m1);
              }  
	       if((m2 <= n)  && (isnan(u[i][m2])))
                  {
                    u[i][m2] = interpH(u, ell, n0, interp_points, i, m2);
                  }       
                  
                }
            }
        } /*terminou a interpolaĆ§Ć£o nas colunas*/

      


    }/*termina o ciclo em j*/
}

/*=======================================================================*/
int local_scalex(double** u,  int k, int m, int ell,  int n0)
/*=======================================================================*/
/* Retorna o minimo entre a primeira casa  ą direita ou ą esquerda de k em que u nao e NAN e o nivel do ponto */

{
 int n = (n0) * (1 << ell);
 int ii; 
 int DistxLeft;
 int DistxRight;
 int esc,level;
 int inf[3];

    	  ii=k-1;
	 while ((ii>0)&&(isnan(u[ii][m]))) 
	   {
	   ii=ii-1; 
	     }
	    DistxLeft=k-ii;
	     
	     ii=k+1;
	  while ((ii<n) && (isnan(u[ii][m]))) 
	   {ii=ii+1;
	    }
	    DistxRight=ii-k;
	  
	     if (DistxLeft< DistxRight)
	   {esc= DistxLeft;
	   }
	   else
	   {esc= DistxRight;
	   }
	   
	   find_level(inf, k, m, ell); 
	   level= 1<<(ell-inf[0]); 
	    /*fprintf(stderr, "level %d\n", level);*/
           if (esc< level) return esc;
	   return level;
	   
	   

   
}

/*=======================================================================*/
int local_scaley(double** u,  int k, int m, int ell, int n0)
/*=======================================================================*/
/* Retorna o minimo entre a primeira casa  abaixo ou acima de m em que u nao e NAN
e o nivel do ponto*/

{
 int n = (n0) * (1 << ell);
 int ii; 
 int DistyDown;
 int DistyUp;
 int esc,level;
 int inf[3];


     
    	 
    	  ii=m-1;
	  while ((ii>0)&&(isnan(u[k][ii]))) 
	   {ii=ii-1; 
	    }
	    DistyDown=m-ii;
	    
    	  ii=m+1;
	  while ((ii<n) && (isnan(u[k][ii]))) 
	   {ii=ii+1; 
	    }
	    DistyUp=ii-m;
	  
	     if (DistyDown < DistyUp)
	   {esc= DistyDown;
	   }
	   else
	   {esc= DistyUp;
	   }
	   
	   find_level(inf, k, m, ell); 
	   level= 1<<(ell-inf[0]); 
	    /*fprintf(stderr, "level %d\n", level);*/
           if (esc< level) return esc;
	   
	   return level;
      
   
}




/*=======================================================================*/
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*/
}

/* UTILITY ROUTINES */

void wait_for_user(void)
  { char buf[10];
    fputs("\nHit return to continue", stderr);
    fgets(buf, sizeof(buf), stdin);
  }

#define BUFSZ 1000

double get_double_from_user(char *prompt, double defval)
  {
    char buf[BUFSZ];
    fprintf(stderr,"Enter %s", prompt);
    if (! isnan(defval)) { fprintf(stderr," [%g]", defval); }
    fprintf(stderr,": ");
    fgets(buf, BUFSZ, stdin);
    char *data = buf;
    while (((*data) == ' ') || ((*data) == '\011') || ((*data) == '\015')) { data++; }
    if (((*data) != '\n') && ((*data) != 0))
      { /* Parse input as integer: */
        char *rest = NULL;
        double dval = strtod(data, &rest); 
        /* Check number syntax: */
        while(((*rest) == ' ') || ((*rest) == '\011') || ((*rest) == '\015')) { rest++; }
        if ((*rest) != '\n')
          { fprintf(stderr, "** invalid character '%c' = '\\%03o'\n", (*rest), (*rest));
            exit(-1);
          }
        return dval;
      }
    else
      { /* Empty input: */
        return defval;
      }
  }

int get_int_from_user(char *prompt, int defval)
  {
    /* All 32-bit ints are exactly representable as {double}s, so: */
    double dval = get_double_from_user(prompt, (double)defval); 
    /* Check whether value is an exact integer: */
    int ival = (int)dval;
    if (dval != ((double)ival)) 
      { fprintf(stderr, "** value '%g' must be an integer\n", dval);
        exit(-1);
      }
    return ival;
  }

void fatal(char* s)
  { fprintf(stderr, "ERROR: %s\n", s);
    perror(s);
    exit(-1);
  }

double** dvector(int nl, int nc)
  { int i;
    double** dv = (double**) malloc(nl*sizeof(double*));
    if (dv == NULL)
      { fprintf(stderr, "ERROR: allocation failure in dvector(%d, %d)\n", nl, nc);
        exit(-1);
      }
    for(i = 0; i < nl; i++)
      { dv[i] = (double*)malloc(nc*sizeof(double));
        if (dv[i] == NULL)
          { fprintf(stderr, "ERROR: allocation failure in dvector(%d, %d)\n", nl, nc);
            exit(-1);
          }
      }
    return dv;
  }
void libera(int nl, int nc, double **A)
  { int i;
    for(i = 0; i < nl; i++)
       { free(A[i]);
        }
    free(A);
  }