/* Last edited on 2024-11-20 06:06:37 by stolfi */ 

/*
 * spr2D.c
 * 29/11/2006
 *
 * Demonstrates sparse decompositions based on iterated 
 * interpolation in 2D   Gaussian
 * Available geometries: empty domain, closed box, open horn.
 */

/* We need to define _GNU_SOURCE to get {asprintf}. */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "int2d.h"
#include "mwplot.h"

/* The maximum number of samples that we can have: */
#define NPx 100
#define NPy 100

#define SIGMA_MAX (1e7)
  /* Nominal conductivity {sigma} of metals. */

void fatal(char* s);

double** dvector(int nl, int nc);
  /* Allocates an array with {nl} rows and {nc} columns. */

void wait_for_user(void);
  /* Writes "Hit return to continue" and waits for user to do so. */

/* NON-HOMOGENEOUS MEDIUM

  Non-homogeneous media are simulated by setting 
  the appropriate elements of the {sigma} field. */
  
void setup_geometry_box(int n, double **sigma);
  /* Draws into {sigma} the outline of a microwave horn, with open
    mouth and impervious sides. */

void setup_geometry_horn(int n, double **sigma);
  /*  Draws into {sigma} a square box, with empty interior
   and impervious walls. */

void show_field
  ( mwplot_t *wp, 
    int nx, 
    int ny, 
    double dx, 
    double dy,
    double **A, 
    double vMax,
    char *zLabel,
    double time
  );
  /* Plots on {wp} the field {A}, assumed to be an array of samples
    with {nx} rows and {ny} columns, with step sizes {dx} and {dy},
    respectively. The vertical scale is set assuming the values of {A}
    are in {[-vMax _ +vmax]}.
    
    The Z axis will be labeled with {zLabel}, and the plot will have a
    title showing the specified {time}, assumed to be in nanoseconds.
    
    If the plotter {wp} is configured to write to disk instead of the
    screen, the plot file will be called "{zLabel}-{NNNNNN}.{ext}",
    where {NNNNNN} is the {time} value rounded to integer, and {ext}
    is the appropriate extension ("eps", "png", etc.) */
  
int main()
  {
    /* Fundamental constants */
    double cc = 2.99792458e8;                 /* speed of light in free space */
    double muz = 4.0 * M_PI * 1.0e-7;         /* permeability of free space */
    double epsz =1.0 / (cc*cc*muz);           /* permittivity of free space */
    /* double aimp = sqrt(muz/epsz);  */      /* wave impedance in free space */
    int interp_points; /* 2 or 4 point interpolation (linear or cubic) */
    double eps;    /* error / threshold for the sparse representation */
    double eps1;   /* error / threshold for the sparse representation */
    int n0;        /* number of samples in the coarse level */
    int ell; 
    int ell_atual; /* actual number of decomp levels with significant wavelet coeffs */
    int level_max; /* number of decomposition levels */
    int n;   
    int npontos;   

    char buf[1024];
    int i, j, k;

    /* Starts the gnuplot process, sets some unchanging plot options: */
    mwplot_t *wp = mwplot_start_gnuplot(0);
    mwplot_set_terminal(wp, /*x11*/ 0);
    mwplot_set_palette(wp, /*blue--red*/ 20);

    n0 = 8;
    printf("Enter the number of samples in the coarse level, up to %d: [%d] ", NPx, n0);
    fgets(buf, sizeof(buf), stdin);
    if (buf[0] != '\n') n0 = atoi(buf);

    ell = 4;
    printf("Enter the number of decompositions levels: [%d] ", ell);
    fgets(buf, sizeof(buf), stdin);
    if (buf[0] != '\n') ell = atoi(buf);
    ell_atual = ell;
    eps = 1e-6;
    printf("Epsilon: [%g] ", eps);
    fgets(buf, sizeof(buf), stdin);
    if (buf[0] != '\n') eps = atof(buf);
    eps1=eps/10;

    interp_points = 4;
    printf("Enter the number of interpolation points: [%d]", interp_points);
    fgets(buf, sizeof(buf), stdin);
    if (buf[0] != '\n') interp_points = atoi(buf);

    n = n0 * (1<<ell);

    double **ex=dvector(n+1, n+1);     /* the function f(x) */
    double **ey=dvector(n+1, n+1);     /* the function f(x) */
    double **hz=dvector(n+1, n+1); 
    double **dyex=dvector(n+1, n+1);   /* sparse form of f(x) */
    double **dxey=dvector(n+1, n+1);   /* sparse form of f(x) */
    double **dyhz=dvector(n+1, n+1);   /* sparse form of f(x) */
    double **dxhz=dvector(n+1, n+1);   /* sparse form of f(x) */
    double **az=dvector(n+1, n+1);     /* hz on full grid */
    double **sigma=dvector(n+1, n+1);  /* geometry */
    /*double **coef=dvector(n+1, n+1);  */         /* coeficiente */

    double dx = 1.0/(n);            /* space increment in x-direction */
    double dy = 1.0/(n);            /* space increment in y-direction */
    double dt = 0.5* dx / cc;       /* time step */
    int nmax =(int)(8.0e-9 / dt);   /* total number of time steps */

    double xpulse, ypulse;   /* Position of initial pulse. */
    double wxpulse, wypulse; /* X-width and Y-width of initial pulse. */

    /* define the medium field {sigma} and the position and shape of initial pulse: */
    int geometry = 2;
    switch(geometry)
      {
        case 0:
          fprintf(stderr, "geometry: homogeneous (vacuum)\n");
          xpulse = 0.500; wxpulse = sqrt(1.0/150.0);
          ypulse = 0.500; wypulse = sqrt(1.0/150.0);
          break;

        case 1:
          fprintf(stderr, "geometry: closed conductive box\n");
          setup_geometry_box(n, sigma);
          xpulse = 0.500; wxpulse = 1.0/50.0;
          ypulse = 0.500; wypulse = 1.0/50.0;

          break;

        case 2:
          fprintf(stderr, "geometry: conductive horn\n");
          setup_geometry_horn(n, sigma);
          xpulse = 0.375; wxpulse = 1.0/50.0;
          ypulse = 0.500; wypulse = 1.0/50.0;
          break;

        default:
          fprintf(stderr, "invalid geometry code\n");
          exit(-1);
      }
    
    /* show the {sigma} field: */
    show_field(wp, n, n, dx, dy, sigma, SIGMA_MAX, "SIGMA", 0.0);
    wait_for_user();

    /* initialize {hz} with a Gaussian pulse, except where {sigma} is nonzero. */
    for (i=0; i <=n; i++)
      for (j=0; j <=n; j++)
        { double vx = (i*dx - xpulse)/wxpulse;
          double vy = (j*dy - ypulse)/wypulse;
          hz[i][j]= (sigma[i][j] == 0 ? exp(-vx*vx)*exp(-vy*vy) : 0.0);
          ex[i][j] = ey[i][j] = 0.0;
          dyex[i][j] = dxey[i][j] = 0.0;
          dxhz[i][j] = dyhz[i][j] = 0.0;
          /*coef[i][j]=(1.0 -(dt*sigma[i][j])/2.0*epsz)/(1.0 +(dt*sigma[i][j])/2.0*epsz);*/
        }

    show_field(wp, n, n, dx, dy, hz, 0.5, "Hz", 0.0);
    wait_for_user();

    spr_data_to_sparseH(hz, ell, n0, interp_points, eps1);
    spr_refineH(hz, ell, n0,  interp_points);

    /* BEGIN TIME-STEPPING LOOP */
    for (k=1; k<= nmax; k++)
      {
        /* Calculo de Hz*/
        spr_data_to_sparseEx(ex, ell, n0, interp_points, eps);
        spr_data_to_sparseEy(ey, ell, n0, interp_points, eps);

        spr_refineEx(ex, ell, n0,  interp_points);
        spr_refineEy(ey, ell, n0,  interp_points);

        /* Makes sure that {ex,ey,hz} are all NAN or all valid: */
        for (i=0; i <=n; i++)
          for (j=0; j <=n; j++)
            { if ((isnan(ex[i][j])) && ((!isnan(ey[i][j]))||(!isnan(hz[i][j]))))
                ex[i][j]=interpEx(ex, ell, n0, interp_points, i, j); 

              if ((isnan(ey[i][j])) && ((!isnan(hz[i][j]))||(!isnan(ex[i][j]))))
                ey[i][j]=interpEy(ey, ell, n0, interp_points, i, j); 

              if ((isnan(hz[i][j])) && ((!isnan(ey[i][j]))||(!isnan(ex[i][j]))))
                hz[i][j]=interpH(hz, ell, n0, interp_points, i, j); 
            }

        /* Comparando o numero de pontos usados*/ 
        npontos=0;
        for (i=0; i <=n; i++)
          for (j=0; j <=n; j++)
            { if (!isnan(hz[i][j]))
                npontos=npontos+1;
            }
        fprintf(stderr,  "sparse grid: %d points\n", npontos);
        npontos=(n+1)*(n+1);
        fprintf(stderr,  "full grid: %d points\n", npontos);
        /*fim da comparacao*/

        derivYEx(ex, ell, n0, interp_points, dx, dyex, hz);
        derivXEy(ey, ell, n0, interp_points, dx, dxey, hz);

        /*fazendo o update do hz*/
        for (i=0; i <=n; i++)
          for (j=0; j <=n; j++)
            { if (!isnan(hz[i][j]))
                { hz[i][j] = hz[i][j]+dt*(dyex[i][j]-dxey[i][j])/muz; }
            }
        spr_data_to_sparseH(hz, ell, n0, interp_points, eps1); 
        spr_refineH(hz, ell, n0,  interp_points);

        /* derivando o campo {hz} */
        derivYH(hz, ell, n0, interp_points, dx, dyhz, ex);
        derivXH(hz, ell, n0, interp_points, dx, dxhz, ey);

        /* fazendo o update do ex */
        for (i=0; i <=n; i++)
          for (j=0; j <=n; j++)
            { if (!isnan(ex[i][j]))
                { double sij = sigma[i][j];
                  double C = dt*sij/(2.0*epsz);
                  ex[i][j] = ((1-C) * ex[i][j] + (dt/epsz) * dyhz[i][j])/(1+C);
                }
            }

        /*fazendo o update do ey*/
        for (i=0; i <=n; i++)
          for (j=0; j <=n; j++)
            { if (!isnan(ey[i][j]))
                { double sij = sigma[i][j];
                  double C = dt*sij/(2.0*epsz);
                  ey[i][j] = ((1-C) * ey[i][j] - (dt/epsz) * dxhz[i][j])/(1+C);
                }
            }

        level_max = find_level_max(ey, ell, n0, ell_atual);
        fprintf(stderr, "level_max = %d\n", level_max);

        /* Copy {hz} to {az} and remove NaNs for plotting: */
        for (i=0; i <= n; i++)
          for (j=0; j <= n; j++)
            az[i][j]= hz[i][j];
        spr_sparse_to_dataH(az, ell, n0, interp_points); 

        /* Plot the magnetic field {hz}: */
        show_field(wp, n, n, dx, dy, az, 0.5, "Hz", k*dt*1e9);
        /* wait_for_user(); */ 

      }/* fim da evolucao no tempo */ 

    mwplot_stop_gnuplot(wp);
    return 0;
  }

void setup_geometry_box(int n, double **sigma)
  {
    double sigma1 = SIGMA_MAX; /* Conductivity of box walls. */

    int nn4 = n/4;
    int i, j;
    
    /* Clear all elements ({malloc} does not do it!): */
    for (i = 0; i < n; i++)
      for (j = 0; j < n; j++)
        { sigma[i][j] = 0; }

    /* Draw the box: */
    for (i = nn4; i <= n-nn4; i++)
      { sigma[nn4][i]= sigma1;
        sigma[n-nn4][i]= sigma1;
        sigma[i][nn4]=sigma1;
        sigma[i][n-nn4]=sigma1;
      }
    for (i= nn4+1; i <= n-nn4-1; i++)
      { sigma[nn4+1][i]= sigma1;
        sigma[n-nn4-1][i]= sigma1;
        sigma[i][nn4+1]=sigma1;
        sigma[i][n-nn4-1]=sigma1;
      }
    for (i = nn4+2; i <= n-nn4-2; i++)
      { sigma[nn4+2][i]= sigma1;
        sigma[n-nn4-2][i]= sigma1;
        sigma[i][nn4+2]=sigma1;
        sigma[i][n-nn4-2]=sigma1;
      }
  }

void setup_geometry_horn(int n, double **sigma)
  {
    double sigma1 = SIGMA_MAX; /* Conductivity of horn walls */

    /* The horn stops at x = 3/4. */
    int nn4 = n/4;
    int nn = 3*n/8;         
    int kn = n/16;   /* no. de degraus de largura w*dx; k=n/4*w, w=4. */

    int i, j, ver, horn, cont;

    /* Clear all elements ({malloc} does not do it!): */
    for (i = 0; i < n; i++)
      for (j = 0; j < n; j++)
        { sigma[i][j] = 0; }

    /* Draw the horn: */
    for (i = nn; i <= n-nn; i++)  /* fundo da corneta */
      { sigma[nn4][i]= sigma1;
        sigma[nn4-1][i]= sigma1;
        sigma[nn4+1][i]= sigma1;
      }
    for (i = nn4-1; i <= n/2; i++)  /* horn parallel sides */
      { sigma[i][nn]= sigma1;    /* bottom side (3 layers) */
        sigma[i][nn-1]= sigma1;
        sigma[i][nn+1]= sigma1;
        sigma[i][n-nn]= sigma1;  /* top side (3 layers) */
        sigma[i][n-nn+1]= sigma1;
        sigma[i][n-nn-1]= sigma1;
      }
    horn=n/2; /* Start of horn's opening */
    ver=nn-1;
    for (cont = 1; cont <= kn; cont++)  /* counting the steps */
      { for (i = horn; i <= horn+4; i++)
          { sigma[i][ver]= sigma1;
            sigma[i][ver-1]= sigma1;
            sigma[i][ver+1]= sigma1;
            sigma[i][n-ver]= sigma1;
            sigma[i][n-ver+1]= sigma1;
            sigma[i][n-ver-1]= sigma1;
          }
        ver=ver-1;  /* descending one step */
        horn=horn+4;            
      }           
  }

void show_field
  ( mwplot_t *wp, 
    int nx, 
    int ny, 
    double dx, 
    double dy,
    double **A, 
    double vMax,
    char *zLabel,
    double time
  )
  { 
    mwplot_set_bivariate_plot_style(wp, /*surface with palette*/ 10);
    mwplot_set_axis_labels(wp, "X", "Y", zLabel);
    mwplot_set_ranges(wp, -0.5*dx, (nx+0.5)*dx, -0.5*dy, (ny+0.5)*dy, -vMax, vMax);
    /* Make up a {title} string: */
    char *title = jsprintf("Time = %g ns", time);
    /* Make up a file name, just in case: */
    char *fileName = jsprintf(&fileName, "%s-%06d", zLabel, (int)floor(time + 0.5));
    mwplot_plot_array(wp, fileName, nx, ny, dx, dy, A, title);
    free(title);
    free(fileName);
  }

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

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(nc*sizeof(double*));
    if (dv == NULL)
      {
        fprintf(stderr, "ERROR: allocation failure in dvector(%d, %d)\n", nl, nc);
        exit(-1);
      }

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