/* última atualização 23/06/2011 */
/* Tenho que ajustar o programa para poder escolher o tempo final e ajustar o número de passos ao tempo final. 

/* 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 <assert.h>
#include <values.h>
#include <time.h>

#include "marktime.h"
#include "Int2Dlocal.h"
#include "mwplot.h"
#include "deriv.h"

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

#define SIGMA_MAX (1.0e9)
  /* Nominal conductivity {sigma} of metals. */
  
void setup_initial_state
  ( int geom, 
    int nx, 
    int ny, 
    double dx, 
    double dy, 
    double **sigma, 
    double **hzx, 
    double **hzy,
    double **ex, 
    double **ey
  );
  /* Sets the conductivity map {sigma} according to the
    code {geom} that describes the domain's geometry.  Currently:
    
      0 -> homogeneous non-conductive medium (vacuum).
      1 -> closed square conductive box.
      2 -> conductive microwave horn.
      
    Also sets the initial state of the magnetic field {hz} and of
    the electric field {ex,ey} to values appropriate to that 
    geometry. 
  */

void setup_geometry_box(int nx, int ny, double **sigma);
  /* Draws into {sigma[0..nx][0..ny]} the outline of a microwave horn,
    with open mouth and highly conductive sides. */

void setup_geometry_horn(int nx, int ny, double **sigma);
  /* Draws into {sigma[0..nx][0..ny]} a closed square box, with empty
    interior and highly conductive walls. */

/* DATA PLOTTING */

void show_field
  ( mwplot_t *wp, 
    int nx, 
    int ny, 
    double dx, 
    double dy,
    double **A, 
    double vMag,
    char *zLabel,
    double time,
    int espera
  );
  /* generates on {wp} a 3D graph of 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 {[-vMag _ +vMag]}.
    
    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.).
    
    If {espera} ir TRUE, waits for user Ok to continue. */
  
  
void show_grid
  ( mwplot_t *wp, 
    int nx, 
    int ny, 
    double dx, 
    double dy,
    double **A, 
    double vMag,
    char *zLabel,
    double time,
    int espera
  );
  /* Generates on {wp} a 2D plot of the matrix {A}, assumed to be an
    array of samples with {nx} rows and {ny} columns, with step sizes
    {dx} and {dy}, respectively. The color scale is set assuming the
    values of {A} are in {[-vMag _ +vMag]}.
    
    The plot will be a grid of colore squares. It 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.).
    
    If {espera} ir TRUE, waits for user Ok to continue. */
    
void print_dadosEx(int nx,int ny, int k, double **A); // prints matrix Ex at time t
void print_dadosEy(int nx,int ny, int k, double **A); // prints matrix Ey at time t
void print_dadosHz(int nx,int ny, int k, double **A); // prints matrix Hz at time t

void get_grid_structure(int nx, int ny, double **F, double **A);
  /* Sets {A[0..nx][0..ny]} to a map showing the occupied elements
    of the sparse arrray {F[0..nx][0..ny]}.  Specifically.
    each {A[i][j]} is set to 0.0 if {F[i][j]} is NAN, and
    to 1 otherwise. The array {A} can then be displayed with
    {show_grid}. */

void get_field_range(int nx, int ny, double **F, double *vLo, double *vHi, double *vMag);
  /* Stores in {*vLo} and {*vHi}, respectively, the minimum and
    maximum values seen in {F[0..nx][0..ny]} Also stores in {*vHi}
    the maximum absolute value among all elements -- that is,
    {max{fabs(vHi),fabs(vLo)}}. */

void print_grid_statistics(FILE *npf, double t, int nx, int ny, double **A);
  /* Counts the number {nvalid} of valid (not NaN) entries in {A[0..nx][0..ny]},
    compares with total number of entries {ntotal=(nx+1)*(ny+1)}.
    Also writes a line into file {npf} with the time {t}
    and the fill ratio {nvalid/ntotal}. */

void main_spr
  ( mwplot_t *wpf, 
    mwplot_t *wpg, 
    FILE *npf,
    int n0,
    int ell,
    int N, 
    int interp_points,
    double dx, 
    double dy,
    double dt,
    int nsteps,
    double **hzx, 
    double **hzy, 
    double **ex, 
    double **ey, 
    double **sigma,
    double **tplt
  );

void main_pml
  ( mwplot_t *wpf, 
    mwplot_t *wpg, 
    FILE *npf,
    int n0,
    int ell,
    int N, 
    int interp_points,
    double dx, 
    double dy,
    double dt,
    int nsteps,
    double **hzx, 
    double **hzy, 
    double **ex, 
    double **ey, 
    double **sigma,
    double **tplt
  );

/* MAIN PROGRAM */

/* Fundamental constants */
#define cc (2.99792458e8)          /* speed of light in free space */
#define muz (4.0 * M_PI * 1.0e-7)   /* permeability of free space */
#define epsz (1.0 / (cc*cc*muz))    /* permittivity of free space */

int main(int argc, char **argv)
  {
   fprintf(stderr, "Max samples in coarsest level = %d\n", NPx);
   int n0 = get_int_from_user("Number of samples in the coarsest level (min 64)", 64);
   int ell = get_int_from_user("Number of decompositions levels (min 2)", 4); 
   int N = n0 * (1<<ell); /* Number of samples in the finest level. */
    
   fprintf(stderr, "n = %d\n", N);

   /* Get order {interp_points} of interpolation (2=linear or 4=cubic): */
   int interp_points = 4;//get_int_from_user("the number of interpolation points", 4);

   double dx = 1.0/N;              /* space increment in x-direction */
   double dy = 1.0/N;              /* space increment in y-direction */
   double dt = 0.7*dx/cc;          /* time step */
   //double dt = 0.35*dx/cc;  
    
   int nsteps = 600;
   double tMax = 1.0;           /* total integration time. */
   fprintf(stderr, "tmax = %.4e  dt = %.4e  nsteps = %d\n", tMax, dt, nsteps);

   double **ex=dvector(N+1, N+1);     /* X component of the elecric field. */
   double **ey=dvector(N+1, N+1);     /* Y component of the electric field. */
   double **hzx=dvector(N+1, N+1);     /* Z component of the magnetic field. */
   double **hzy=dvector(N+1, N+1);     /* Z component of the magnetic field. */
    
   double **tplt=dvector(N+1, N+1);   /* Temporary array for field plots. */
   double **sigma=dvector(N+1, N+1);  /* Conductivity of medium. */

   /* Starts the gnuplot process for 3D field plots: */
   mwplot_t *wpf = mwplot_start_gnuplot(640,480);
   mwplot_set_terminal(wpf,  0); /* write to terminal*/
   //mwplot_set_terminal(wpf, 11);   /* write to disk (eps color).*/ 

   /* Starts the gnuplot process for 2D grid plots */
   mwplot_t *wpg = mwplot_start_gnuplot(640,480); /* 0,0 for debug output */
   //mwplot_set_terminal(wpg,  0); /*write to terminal*/
   mwplot_set_terminal(wpg,  11); /*write to disk*/
   mwplot_set_gridmap_palette(wpg, /*yellow--red*/ 0);

   /* Open file {npf} with plot of matrix fill ratio over time: */
   char *npfname = NULL;
   asprintf(&npfname, "numberOfNAN-%03d.dat", ell);
   FILE *npf = fopen(npfname, "w");
   assert(npf != NULL);

   /* Define the conductivity field {sigma} and the initial pulse: */
   int geometry = get_int_from_user("the conductivity map (domain geometry) code", 0);
   setup_initial_state(geometry, N, N, dx, dy, sigma, hzx, hzy, ex, ey);
    
   /* Show the {sigma} field: */
   //show_grid(wpg, N, N, dx, dy, sigma, SIGMA_MAX, "SIGMA", 0.0, 1);

   int pml = get_int_from_user("0 for SPR, 1 for PML", 1);
   if (pml)
     {
      main_pml(wpf, wpg, npf, n0, ell, N, interp_points, dx, dy, dt, nsteps, hzx, hzy, ex, ey, sigma, tplt);
     }
   else
     {
      main_spr(wpf, wpg, npf, n0, ell, N, interp_points, dx, dy, dt, nsteps, hzx, hzy, ex, ey, sigma, tplt);
     }

   mwplot_stop_gnuplot(wpf);
   mwplot_stop_gnuplot(wpg);
   fclose(npf);
   return 0;
  }

void main_spr
  ( mwplot_t *wpf, 
    mwplot_t *wpg,
    FILE *npf,
    int n0,
    int ell,
    int N, 
    int interp_points,
    double dx, 
    double dy,
    double dt,
    int nsteps,
    double **hzx, 
    double **hzy, 
    double **ex, 
    double **ey,
    double **sigma,
    double **tplt
  )
  { /* double aimp = sqrt(muz/epsz);  */      /* wave impedance in free space */

   int i, j, k;
   
   double eps = get_double_from_user("the error/threshold for electical field sparse representation", 0.01);
   double eps1 = eps/100;  /* error/threshold for the sparse representation. */
   double **dA1=dvector(N+1, N+1);   /* dex/dy - sparse form. */
   double **dA2=dvector(N+1, N+1);   /* dey/dx - sparse form. */
     
   /*===============================================================*/

   /* 2^{-n0} - related to the number of cells/thickness of PML*/
   //int maxNumberOfPMLCells = 64;	
   //maxNumberOfPMLCells;  // To be change 06/11/2009
   int L = N/8;// (1<<ell);  /* PML border width (in steps; actual width is {L*dx}). */
   double d = ((double)L)/((double)N);

   /* Constants: */
   double lnR0 = log(1.0e-7);                /* ln(R(0)) */
   double or = 2.00;   /* order of PML*/
   double eta = sqrt(muz/epsz); /* wave impedance in free space ???*/
   double sigmax = -(or + 1.0)*lnR0/(2.0*eta*d); /* max {sigma} at the PML bdry. */
   double temp_proc;
   double a1x, a1y, a2x, a2y;

   /*definindo as variaveis da regiao de PML*/
   double sig[N+1];
   double sigs[N+1];
   double c1[N+1];
   double c2[N+1];
   double c1s[N+1];
   double c2s[N+1];
   for(i = 0; i < N; i++)
     sig[i] = sigs[i] = c1[i] = c2[i] = c1s[i] = c2s[i] =0.0;
		
		
   /* Set up the PML material coefficients */
   for(i = 0; i < L; i++)
     {	
      sig[i] = sigmax*pow(((L-i)*1.0)/(L),or);
      sigs[i] = sig[i]*muz/epsz;
      c1[i]=exp(-sig[i]*dt/epsz);
      c1[N-i]=c1[i];
      c2[i]=(1.0-c1[i])/sig[i];
      c2[N-i]=c2[i];
      c1s[i]=exp(-sigs[i]*dt/muz);
      c1s[N-i]=c1s[i];
      c2s[i]=(1.0-c1s[i])/sigs[i];
      c2s[N-i]=c2s[i];
     }

   for (i=L; i <= N-L; i++)
    { 
     c1[i] = c1s[i]=1.0;
     c2[i] = dt/epsz;
     c2s[i]= dt/muz;
    }

   /*===============================================================*/
    
   /*double **coef=dvector(N+1, N+1);  */         /* coeficiente */
   /* Find the maximum {hz} value: */
   double hzLo, hzHi, hzMag;
   for(i = 0; i <= N; i++)
    for (j = 0; j <= N; j++)
     tplt[i][j]= hzx[i][j]+hzy[i][j];
   get_field_range(N, N, tplt, &hzLo, &hzHi, &hzMag);
    
   /* Show the initial {hz} field, allow user to choose the style: */
   int pal = 0; /* User-selected palette code. */
   int sty = 0; /* User-selected plot style code. */
   while (1)
     { /* Set the selected palette and style: */
      fprintf(stderr, "using 3D graph palette %d\n", pal);
      mwplot_set_3dgraph_palette(wpf, pal);
      fprintf(stderr, "using 3D graph style %d\n", sty);
      /* Display the field: */
      mwplot_set_3dgraph_style(wpf, sty);
      show_field(wpf, N, N, dx, dy, tplt, hzMag, "Hz", 0.0, 0);
      //print_dadosHz(N, N, 0.0, tplt); // salva condicao inicial em arquivo
      /* Ask user for a new palette and style: */
      pal = get_int_from_user("a new palette code, or -1 to go on", -1);
      if (pal == -1) { break; }
      sty = get_int_from_user("a new style code", sty);
     } 
    
   /* Convert the initial field {hz} to sparse representation: */
   spr_data_to_sparseH(hzx, ell, n0, interp_points, eps1);
   spr_refineH(hzx, ell, n0,  interp_points);
   spr_data_to_sparseH(hzy, ell, n0, interp_points, eps1);
   spr_refineH(hzy, ell, n0,  interp_points);
    
   /* Show the initial sparse grid for {hz}, allow user to choose style: */
   pal = 0;
   sty = 0;
   while (1)
     { /* Set the selected palette and style: */
      fprintf(stderr, "using grid map palette %d\n", pal);
      mwplot_set_gridmap_palette(wpg, pal);
      fprintf(stderr, "using grid map style %d\n", sty);
      /* Display the field: */
      mwplot_set_gridmap_style(wpg, sty);
      get_grid_structure(N, N, hzx, tplt);
      //show_grid(wpg, N, N, dx, dy, tplt, 1.0, "GRID", 0.0, 0);
      /* Ask user for a new palette and style: */
      pal = get_int_from_user("a new palette code, or -1 to go on", -1);
      if (pal == -1) { break; }
      sty = get_int_from_user("a new style code", sty);
     }
   
   /* Clear out the field derivatives {dyex,dxey,dxhz,dyhz}: */
   for (i=0; i <= N; i++)
     for (j=0; j <= N; j++)
       {       
        dA1[i][j] = dA2[i][j] = 0.0;
       }

   double time, tMax=2.0;
   
 
   /* BEGIN TIME-STEPPING LOOP */
   int freq =(1<< ell);
   freq=60;
   MarkTime (NULL);
   printf("começou\n");
   for (k = 1; k <= nsteps; k++)
     {	
      int show = (k % freq == 0) || (k == 1);
      time = k * dt * 1.0e9;
      //printf("o valor de k é %d\n",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(hzx[i][j])) || (!isnan(hzy[i][j])))
             {
	      if (isnan(ex[i][j]))  ex[i][j]=interpEx(ex, ell, n0, interp_points, i, j); 
              if (isnan(ey[i][j]))  ey[i][j]=interpEy(ey, ell, n0, interp_points, i, j); 
              if (isnan(hzx[i][j])) hzx[i][j]=interpH(hzx, ell, n0, interp_points, i, j); 
	      if (isnan(hzy[i][j])) hzy[i][j]=interpH(hzy, ell, n0, interp_points, i, j);  
	     }    
          }

      derivYEx(ex, ell, n0, interp_points, dy, dA2, hzy);
      derivXEy(ey, ell, n0, interp_points, dx, dA1, hzx);
      /* Update the {hz} field: */
      // Free space
      a1x= a1y=1.0;
      a2x= a2y=dt/muz;
      for (i = L; i < N-L; i++)
        for (j = L; j < N-L; j++)
          {
           if (!isnan(hzx[i][j]))
             { 	       
              hzx[i][j]=a1x*hzx[i][j]-a2x*dA1[i][j];
	      hzy[i][j]=a1y*hzy[i][j]+a2y*dA2[i][j]; 
	     }
	  }
      //front     
      for (i = 0; i <= N; i++)
        {
         a1x=c1s[i];
	 a2x=c2s[i];
         for (j = 0; j < L; j++)
           {
            if (!isnan(hzx[i][j]))
	     {
              a1y=c1s[j];
	      a2y=c2s[j];
              hzx[i][j]=a1x*hzx[i][j]-a2x*dA1[i][j];
	      hzy[i][j]=a1y*hzy[i][j]+a2y*dA2[i][j]; 
	     }
	   } 
	} 
      //top
      for (i = 0; i <= N; i++)
	{
	 a1x=c1s[i];
	 a2x=c2s[i];
         for (j = N-L; j <=N; j++)
           {
            if (!isnan(hzx[i][j]))
	     {
              a1y=c1s[j];
	      a2y=c2s[j];
              hzx[i][j]=a1x*hzx[i][j]-a2x*dA1[i][j];
	      hzy[i][j]=a1y*hzy[i][j]+a2y*dA2[i][j]; 
	     }
	   }
	}     
      //left      
      for (i = 0; i < L; i++)
	{
	 a1x=c1s[i];
	 a2x=c2s[i];
         for (j = L; j < N-L; j++)
           {
            if (!isnan(hzx[i][j]))
	      {
               a1y=c1s[j];
	       a2y=c2s[j];
               hzx[i][j]=a1x*hzx[i][j]-a2x*dA1[i][j];
	       hzy[i][j]=a1y*hzy[i][j]+a2y*dA2[i][j]; 
	      }
	   } 
	}
      //right    
      for (i = N-L; i <= N; i++)
	{
         a1x=c1s[i];
	 a2x=c2s[i];
         for (j = L; j < N-L; j++)
           {
            if (!isnan(hzx[i][j]))
	      {
               a1y=c1s[j];
	       a2y=c2s[j];
               hzx[i][j]=a1x*hzx[i][j]-a2x*dA1[i][j];
	       hzy[i][j]=a1y*hzy[i][j]+a2y*dA2[i][j]; 
	      }
	   }  
	}
      spr_data_to_sparseH(hzx, ell, n0, interp_points, eps1); 
      spr_data_to_sparseH(hzy, ell, n0, interp_points, eps1); 
	
      //if (k==nsteps)
      //  { 
      //   get_grid_structure(N, N, hzx, tplt);
      //   show_grid(wpg, N, N, dx, dy, tplt, 1.0, "GRID", time, 0);
      //  }

      spr_refineH(hzx, ell, n0,  interp_points);
      spr_refineH(hzy, ell, n0,  interp_points);

      /* Differentiate the {hzx and hzy} fields: */
      derivYH(hzx, ell, n0, interp_points, dy, dA1, ex);
      derivYH(hzy, ell, n0, interp_points, dy, dA2, ex);
	 
      /* Updating the {ex} field: */
      //free space
      for (i = L; i < N-L; i++)
        for (j = L; j < N-L; j++)
	  {
           if (!isnan(ex[i][j]))
	     {
              a1x=1.0;
	      a2x=dt/epsz;
	      if(sigma[i][j]!=0.0)
	        { 
	         a1x=exp(-sigma[i][j]*dt/epsz);
                 a2x= (1-exp(-sigma[i][j]*dt/epsz))/sigma[i][j];
                }
	      ex[i][j]=a1x*ex[i][j]+a2x*(dA1[i][j]+dA2[i][j]);
	     }	
          }
      // front    
      for (i = 0; i <= N; i++)
        for (j = 0; j < L; j++)
          {
           if (!isnan(ex[i][j]))
             {
              a1x=c1[j];
	      a2x= c2[j];
	      ex[i][j]=a1x*ex[i][j]+a2x*(dA1[i][j]+dA2[i][j]);
	     }
          }
      //top  
      for (i = 0; i <= N; i++)
        for (j = N-L; j <= N; j++)
          {
           if (!isnan(ex[i][j]))
             {
              a1x=c1[j];
	      a2x= c2[j];
	      ex[i][j]=a1x*ex[i][j]+a2x*(dA1[i][j]+dA2[i][j]);
	     }
          }
      //left
      for (i = 0; i < L; i++)
        for (j = L; j < N-L; j++)
          {
           if (!isnan(ex[i][j]))
             {
              a1x=c1[j];
	      a2x= c2[j];
	      ex[i][j]=a1x*ex[i][j]+a2x*(dA1[i][j]+dA2[i][j]);
	     }
          }
      //right  
      for (i = N-L; i <= N; i++)
        for (j = L; j < N-L; j++)
          {
           if (!isnan(ex[i][j]))
             {
              a1x=c1[j];
	      a2x= c2[j];
	      ex[i][j]=a1x*ex[i][j]+a2x*(dA1[i][j]+dA2[i][j]);
	     }	
          }
         
      /* Updating the {ey} field: */
         
      derivXH(hzx, ell, n0, interp_points, dx, dA1, ey);
      derivXH(hzy, ell, n0, interp_points, dx, dA2, ey);
	
      //free
      for (i = L; i < N-L; i++)
        for (j = L; j < N-L; j++)
          {
           if (!isnan(ey[i][j]))
	     {
              a1y=1.0;
	      a2y=dt/epsz;
	      if(sigma[i][j]!=0.0)
	        {
                 a1y=exp(-sigma[i][j]*dt/epsz);
	         a2y= (1-exp(-sigma[i][j]*dt/epsz))/sigma[i][j];
                }
              ey[i][j]=a1y*ey[i][j]-a2y*(dA1[i][j]+dA2[i][j]);
	     }
          }
      //front
      for (i = 0; i <= N; i++)
	{
         a1y=c1[i];
	 a2y= c2[i];
         for (j = 0; j < L; j++)
           {
            if (!isnan(ey[i][j]))
              ey[i][j]=a1y*ey[i][j]-a2y*(dA1[i][j]+dA2[i][j]);
	   }
        }
      // top
      for (i = 0; i <= N; i++)
	{
         a1y=c1[i];
	 a2y= c2[i];
         for (j = N-L; j <= N; j++)
           {
            if (!isnan(ey[i][j]))
              ey[i][j]=a1y*ey[i][j]-a2y*(dA1[i][j]+dA2[i][j]);
	   }			
        }
      //left
      for (i = 0; i < L; i++)
	{
         a1y=c1[i];
	 a2y= c2[i];
         for (j = L; j < N-L; j++)
           {
            if (!isnan(ey[i][j]))
              ey[i][j]=a1y*ey[i][j]-a2y*(dA1[i][j]+dA2[i][j]);
	   }
        }
      //right
      for (i = N-L; i <= N; i++)
	{
         a1y=c1[i];
	 a2y= c2[i];
         for (j = L; j < N-L; j++)
           {
            if (!isnan(ey[i][j]))
              ey[i][j]=a1y*ey[i][j]-a2y*(dA1[i][j]+dA2[i][j]);
	   }
        }

      /* Find max level used in this iteration: */
      // int ell_atual = ell; /* actual num of levels with signif wavelet coeffs. ????*/
      // int level_max = find_level_max(ey, ell, n0, ell_atual);
      	   
      if (show)
	{ 
         //fprintf(stderr, "level_max = %d\n", level_max);
	 /* Plot the magnetic field {hz}: */
         for(i = 0; i <= N; i++)
           for (j = 0; j <= N; j++)
             tplt[i][j]= ex[i][j];
         spr_sparse_to_dataEx(tplt, ell, n0, interp_points); 
         //show_field(wpf, N, N, dx, dy, tplt, NAN, "Ex", time, 0);
         //print_dadosEx(N, N, k, tplt);
         for(i = 0; i <= N; i++)
           for (j = 0; j <= N; j++)
             tplt[i][j]= ey[i][j];
         spr_sparse_to_dataEy(tplt, ell, n0, interp_points); 
         //show_field(wpf, N, N, dx, dy, tplt, NAN, "Ey", time, 0);
         //print_dadosEy(N, N, k, tplt);
         for(i = 0; i <= N; i++)
           for (j = 0; j <= N; j++)
             tplt[i][j]= hzx[i][j]+hzy[i][j];
         spr_sparse_to_dataH(tplt, ell, n0, interp_points); 
	 show_field(wpf, N, N, dx, dy, tplt, NAN, "Hz", time, 0);
         //print_dadosHz(N, N, k, tplt);
        }
            
     }/* fim da evolucao no tempo */ 
      
   MarkTime (&temp_proc);
   printf("Terminou\n");
   ShowTime( temp_proc, "\n tempo gasto: ");
  }

void main_pml
  ( mwplot_t *wpf, 
    mwplot_t *wpg, 
    FILE *npf,
    int n0,
    int ell,
    int N, 
    int interp_points,
    double dx, 
    double dy,
    double dt,
    int nsteps,
    double **hzx, 
    double **hzy, 
    double **ex, 
    double **ey, 
    double **sigma,
    double **tplt
  )
  {
   

    /* dados do PML*/
        int L = N/8;  /* PML border width (in steps; actual width is {L*dx}). */
        double d = ((double)L)/((double)N);

    /* Constants: */
    double lnR0 = log(1.0e-7);                /* ln(R(0)) */
    double or = 2.00;   /* order of PML*/
    double eta = sqrt(muz/epsz); /* wave impedance in free space ???*/
    double sigmax = -(or + 1.0)*lnR0/(2.0*eta*d); /* max {sigma} at the PML bdry. */
  

    /*definindo as variaveis da regiao de PML*/
    double sig[N+1];
    double sigs[N+1];
    
    
    double c1[N+1];
    double c2[N+1];
    double c1s[N+1];
    double c2s[N+1];

    double a1x, a1y, a2x, a2y;
    int i, j, k;
    
    double temp_proc;

    
 

	for(i = 0; i < N; i++)
		sig[i] = sigs[i] = c1[i] = c2[i] = c1s[i] = c2s[i] =0.0;

    double dyhzx_T;
    double dxhzx_T;
    double dyhzy_T;
    double dxhzy_T;
    double dyex_T;
    double dxey_T;

      
    
    for(i = 0; i < L; i++)
      {	
		sig[i] = sigmax*pow(((L-i)*1.0)/(L),or);
		sigs[i] = sig[i]*muz/epsz;
		
		c1[i]=exp(-sig[i]*dt/epsz);
                c1[N-i]=c1[i];
                c2[i]=(1.0-c1[i])/sig[i];
		c2[N-i]=c2[i];
		c1s[i]=exp(-sigs[i]*dt/muz);
                c1s[N-i]=c1s[i];
		c2s[i]=(1.0-c1s[i])/sigs[i];
		c2s[N-i]=c2s[i];
      }

	for (i=L; i <= N-L; i++)
		{ 
			c1[i] = c1s[i]=1.0;
			c2[i] = dt/epsz;
			c2s[i]= dt/muz;
		}

   
   double time;

   /* Graficando a condição inicial do campo magnetico */
   // for(i = 0; i <= N; i++)
   //   for (j = 0; j <= N; j++)
   //     tplt[i][j]= hzx[i][j]+hzy[i][j];}
   // show_field(wpf, N, N, dx, dy, tplt, NAN, "hz", time, 1);

    
    /* BEGIN TIME-STEPPING LOOP */
    int debug = 0;
    int freq =(1<< ell);

    freq=60;    

    MarkTime (NULL);
    printf("comecou\n");
    for (k = 1; k <= nsteps; k++)
      {
        time = k * dt * 1e9;
	/* int esp = (time > 1.2); */
        int show = (k % freq == 0) || (k == 1);

        /* Update {hzx e hzy}: */
        // Free space
        a1x= a1y=1.0;
	a2x= a2y=dt/muz;
        for (i = L; i < N-L; i++)
          for (j = L; j < N-L; j++)
          {  dxey_T= DerivAS(ey, dx, i, j, N, 1);
	     dyex_T= DerivAS(ex, dy, i, j, N, 0); 
             hzx[i][j]=a1x*hzx[i][j]-a2x*dxey_T;
	     hzy[i][j]=a1y*hzy[i][j]+a2y*dyex_T;
          }
        // front
	for (i = 0; i <= N; i++)
        {  a1x=c1s[i];
	   a2x=c2s[i];
	   for (j = 0; j < L; j++)
           {  dxey_T= DerivAS(ey, dx, i, j, N, 1);
	      dyex_T= DerivAS(ex, dy, i, j, N, 0); 
              a1y=c1s[j];
	      a2y=c2s[j];
	      hzx[i][j]=a1x*hzx[i][j]-a2x*dxey_T;
	      hzy[i][j]=a1y*hzy[i][j]+a2y*dyex_T;
           }
        }
        // top
	for (i = 0; i <= N; i++)
	{  a1x=c1s[i];
	   a2x=c2s[i];
           for (j = N-L; j <= N; j++)
           {  dxey_T= DerivAS(ey, dx, i, j, N, 1);
	      dyex_T= DerivAS(ex, dy, i, j, N, 0); 
              a1y=c1s[j];
	      a2y=c2s[j];
	      hzx[i][j]=a1x*hzx[i][j]-a2x*dxey_T;
	      hzy[i][j]=a1y*hzy[i][j]+a2y*dyex_T;
           }
	}
        // left
	for (i = 0; i < L; i++)
	{  a1x=c1s[i];
	   a2x=c2s[i];
           for (j = L; j < N-L; j++)
           {  dxey_T= DerivAS(ey, dx, i, j, N, 1);
	      dyex_T= DerivAS(ex, dy, i, j, N, 0); 
              a1y=c1s[j];
	      a2y=c2s[j];
	      hzx[i][j]=a1x*hzx[i][j]-a2x*dxey_T;
	      hzy[i][j]=a1y*hzy[i][j]+a2y*dyex_T;
           }
	}
        //right
	for (i = N-L; i <= N; i++)
	{  a1x=c1s[i];
	   a2x=c2s[i];
           for (j = L; j < N-L; j++)
           {  dxey_T= DerivAS(ey, dx, i, j, N, 1);
	      dyex_T= DerivAS(ex, dy, i, j, N, 0); 
              a1y=c1s[j];
	      a2y=c2s[j];
	      hzx[i][j]=a1x*hzx[i][j]-a2x*dxey_T;
	      hzy[i][j]=a1y*hzy[i][j]+a2y*dyex_T;
           }
	}

        /* Updating the {ex} e {ey}field: */
        //free space
      	for (i = L; i < N-L; i++)
          for (j = L; j < N-L; j++)
          {  a1x=1.0;
	     a2x=dt/epsz;
	     if(sigma[i][j]!=0.0)
	     { 
	        a1x=exp(-sigma[i][j]*dt/epsz);
	        a2x= (1-exp(-sigma[i][j]*dt/epsz))/sigma[i][j];
             }
             dxhzx_T= DerivS(hzx, dx, i, j, N, 1);
	     dyhzx_T= DerivS(hzx, dy, i, j, N, 0); 
	     dxhzy_T= DerivS(hzy, dx, i, j, N, 1);
	     dyhzy_T= DerivS(hzy, dy, i, j, N, 0); 
             ex[i][j]=a1x*ex[i][j]+a2x*(dyhzx_T+dyhzy_T);
             ey[i][j]=a1x*ey[i][j]-a2x*(dxhzx_T+dxhzy_T);
          }
        //front
	for (i = 0; i <= N; i++)
        {  a1y=c1[i];
	   a2y=c2[i];
	   for (j = 0; j < L; j++)
           {  a1x=c1[j];
	      a2x=c2[j];
	      dxhzx_T= DerivS(hzx, dx, i, j, N, 1);
	      dyhzx_T= DerivS(hzx, dy, i, j, N, 0); 
	      dxhzy_T= DerivS(hzy, dx, i, j, N, 1);
	      dyhzy_T= DerivS(hzy, dy, i, j, N, 0); 
              ex[i][j]=a1x*ex[i][j]+a2x*(dyhzx_T+dyhzy_T);
              ey[i][j]=a1y*ey[i][j]-a2y*(dxhzx_T+dxhzy_T);
           }
        }
        //top
	for (i = 0; i <= N; i++)
        {  a1y=c1[i];
	   a2y=c2[i];
	   for (j = N-L; j <= N; j++)
           {  a1x=c1[j];
	      a2x= c2[j];
	      dxhzx_T= DerivS(hzx, dx, i, j, N, 1);
	      dyhzx_T= DerivS(hzx, dy, i, j, N, 0); 
	      dxhzy_T= DerivS(hzy, dx, i, j, N, 1);
	      dyhzy_T= DerivS(hzy, dy, i, j, N, 0); 
              ex[i][j]=a1x*ex[i][j]+a2x*(dyhzx_T+dyhzy_T);
              ey[i][j]=a1y*ey[i][j]-a2y*(dxhzx_T+dxhzy_T);
           }
        }
        //left
	for (i = 0; i < L; i++)
        {  a1y=c1[i];
	   a2y=c2[i];
	   for (j = L; j < N-L; j++)
           {  a1x=c1[j];
	      a2x= c2[j];
	      dxhzx_T= DerivS(hzx, dx, i, j, N, 1);
	      dyhzx_T= DerivS(hzx, dy, i, j, N, 0); 
	      dxhzy_T= DerivS(hzy, dx, i, j, N, 1);
	      dyhzy_T= DerivS(hzy, dy, i, j, N, 0); 
              ex[i][j]=a1x*ex[i][j]+a2x*(dyhzx_T+dyhzy_T);
              ey[i][j]=a1y*ey[i][j]-a2y*(dxhzx_T+dxhzy_T);
           }
        }
        //right
	for (i = N-L; i <=N; i++)
        {  a1y=c1[i];
	   a2y=c2[i];
	   for (j = L; j < N-L; j++)
           {  a1x=c1[j];
	      a2x= c2[j];
	      dxhzx_T= DerivS(hzx, dx, i, j, N, 1);
	      dyhzx_T= DerivS(hzx, dy, i, j, N, 0); 
	      dxhzy_T= DerivS(hzy, dx, i, j, N,  1);
	      dyhzy_T= DerivS(hzy, dy, i, j, N,  0); 
              ex[i][j]=a1x*ex[i][j]+a2x*(dyhzx_T+dyhzy_T);
              ey[i][j]=a1y*ey[i][j]-a2y*(dxhzx_T+dxhzy_T);
           }
        }
        /* Mostrar grafico ou gravar arquivo */
         if (show)
         { 
          for(i = 0; i <= N; i++)
           for (j = 0; j <= N; j++)
	     tplt[i][j]= hzx[i][j]+hzy[i][j];
	  show_field(wpf, N, N, dx, dy, tplt, NAN, "hz", time, 0); 
          //print_dadosHz(N,N, k, tplt);
          //print_dadosEx(N,N, k, ex);
          //print_dadosEy(N,N, k, ey);
         }
     
   } /* fim da evolucao no tempo */
     	 
 	
     MarkTime (&temp_proc);
     printf("terminou\n");
     ShowTime( temp_proc, "\n tempo gasto: ");
 }
   
/* IMPLEMENTATIONS */
    
void setup_initial_state
  ( int geom, 
    int nx, 
    int ny, 
    double dx, 
    double dy, 
    double **sigma, 
    double **hzx, 
    double **hzy, 
    double **ex, 
    double **ey
  )
  {
    double xpulse, ypulse;   /* Position of initial pulse. */
    double wxpulse, wypulse; /* X-width and Y-width of initial pulse. */

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

    switch(geom)
      {
        case 1:
          fprintf(stderr, "geometry: closed conductive box\n");
          setup_geometry_box(nx, ny, sigma);
          xpulse = 0.328125; wxpulse = 1.0/50.0;
          ypulse = 0.500000; wypulse = 1.0/50.0;
          //xpulse = (84.0/256.0); wxpulse = 1.0/200.0;
          //ypulse = 0.500; wypulse = 1.0/200.0;
          break;

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

        default:
          fprintf(stderr, "invalid geometry code %d, assuming 0\n", geom);
          geom = 0;

        case 0:
          fprintf(stderr, "geometry: homogeneous non-conductive\n");
          xpulse = 0.500; wxpulse = sqrt(1.0/3500.0);
          ypulse = 0.500; wypulse = sqrt(1.0/3500.0);
          break;
      }

    /* initialize {hz} with a Gaussian pulse, except where {sigma} is nonzero. */
    for (i = 0; i <= nx; i++)
      for (j = 0; j <= ny; j++)
        {ex[i][j] = ey[i][j] = 0.0;
         hzx[i][j]= 0.0;
	 hzy[i][j]= 0.0;
         double vx = (i*dx - xpulse)/wxpulse;
         double vy = (j*dy - ypulse)/wypulse;
	 hzx[i][j]= 0.5*exp(-vx*vx)*exp(-vy*vy);
	 hzy[i][j]= 0.5*exp(-vx*vx)*exp(-vy*vy);
         }
    
    //for (i = 16; i <= 113; i++)  
    //  for (j = 16; j <= 113; j++) 
    //    { double vx = (i*dx - xpulse)/wxpulse;
    //      double vy = (j*dy - ypulse)/wypulse;
    //  hzx[i][j]= 0.5*exp(-vx*vx)*exp(-vy*vy);
    //  hzy[i][j]= 0.5*exp(-vx*vx)*exp(-vy*vy);
    //    }
  }
  
void setup_geometry_box(int nx, int ny, double **sigma)
  { 

    /* The following code assumes a square grid, so: */
    assert(nx == ny);
    int N = nx;
    
    int nn1 = 15*N/64;
    int nn2 = 23*N/64;
    int nn3 = 29*N/64;
    int nn4 = 35*N/64;
    int nn5 = 37*N/64;  
    double dx = 1.0/N;   
    double ix;  
    double jy;   
    int i, j;
    
    /* Draw the box: */
    double sigma1 = SIGMA_MAX; /* Conductivity of box walls. */
    
     /* Draw the box: */
    for (j = nn2; j <= N-nn2; j++)  /* Parede esquerda */
      for (i = 0; i <= N/2; i++) 
      { ix=i*dx;
        if (fabs(ix-0.25) <=0.015625)
          sigma[i][j]= sigma1;
      }

    /* Paredes interiores*/
    for (j = nn2; j <= nn3; j++)  /* Parte de baixo da Primeira parede*/
      for (i = nn1; i <= N/2; i++) 
        { ix=i*dx;
          if (fabs(ix-0.40625) <=0.015625)
            sigma[i][j]= sigma1;
        }
      
    for (j = nn4; j <= N-nn2; j++)  /* Parte de cima da Primeira parede */
      for (i = nn1; i <= N/2; i++) 
      { ix=i*dx;
        if (fabs(ix-0.40625) <=0.015625)
          sigma[i][j]= sigma1;
       }
      
    for (j = nn2; j <= N/2; j++)  /* Parte de baixo da Segunda parede*/
      for (i = N/2; i <= N-nn1; i++) 
      { ix=i*dx;
        if (fabs(ix-0.59375) <=0.015625)
          sigma[i][j]= sigma1;
      }
      
    for (j = nn5; j <= N-nn2; j++)  /* Parte de cima da Segunda parede */
      for (i = N/2; i <= N-nn1; i++) 
      { ix=i*dx;
        if (fabs(ix-0.59375) <=0.015625)
          sigma[i][j]= sigma1;
       }
      
    /* Paredes horizontais. */ 
    for (i = nn1; i <= N-nn1; i++)  
      for (j = 0; j <= N/2; j++)  
      {  jy=j*dx;
         if (fabs(jy-0.375) <=0.015625)
          { sigma[i][j]= sigma1;  
	    sigma[i][N-j]= sigma1; 
	  }             
      }    
 }

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

    /* The following code assumes a square grid, so: */
    assert(nx == ny);
    int N = nx;
    double dx = 1.0/N;   
    double ix;  
    double jy;   
    double yx; 
    
    /* The horn stops at x = 3/4. */
    int nn4 = 15*N/64;
    int nn = 23*N/64;         
   
    int i, j;

    /* Draw the horn: */
    for (j = nn; j <= N-nn; j++)  /* Bottom of horn */
      { for (i = 0; i <= N/2; i++) 
        { ix=i*dx;
          if (fabs(ix-0.25) <=0.015625)
            { sigma[i][j]= sigma1;
            }
         }
       }	 
     
     for (i = nn4; i <= N/2; i++)  /* Horn's parallel sides. */
       {
        for (j = 0; j <= N/2; j++)  
            {  jy=j*dx;
               if (fabs(jy-0.375) <=0.015625)
                { sigma[i][j]= sigma1;  
	          sigma[i][N-j]= sigma1; 
	        }             
	    }    
        }
        
       for (i = N/2; i <= 3*N/4; i++) /* Start of horn's opening */
	   { ix=i*dx;
	      yx= -(ix-0.5)/2.0 + 0.375;
	      for (j = 0; j <= N/2; j++)  
                { jy=j*dx;
                  if (fabs(jy-yx) <=0.015625 )
	            {  sigma[i][j]= sigma1; 
		       sigma[i][N-j]= sigma1; 
		    }
	        }
	     }   	    
        }
  
 
void show_field
  ( mwplot_t *wp, 
    int nx, 
    int ny, 
    double dx, 
    double dy,
    double **A, 
    double vMag,
    char *zLabel,
    double time,
    int espera
  )
  { mwplot_set_3dgraph_style(wp, /*surface with palette*/ 0);
    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, -vMag, vMag);
    /* Make up a {title} string: */
    char *title = NULL;
    asprintf(&title, "Time = %g ns", time);
    /* Make up a file name, just in case: */
    char *fileName = NULL;
    asprintf(&fileName, "%s-%06d", zLabel, (int)floor(1000*time + 0.5));
    mwplot_plot_3dgraph(wp, fileName, nx, ny, dx, dy, A, title);
    if (espera) { wait_for_user(); }
    free(title);
    free(fileName);
  }

void show_grid
  ( mwplot_t *wp, 
    int nx, 
    int ny, 
    double dx, 
    double dy,
    double **A, 
    double vMag,
    char *zLabel,
    double time,
    int espera
  )
  { 
   mwplot_set_gridmap_style(wp, /*default*/ 0);
    mwplot_set_axis_labels(wp, "X", "Y", NULL);
    mwplot_set_ranges(wp, 0, nx+1, 0, ny+1, -vMag, vMag);
    /* Make up a {title} string: */
    char *title = NULL;
    asprintf(&title, "Time = %g ns", time);
    /* Make up a file name, just in case: */
    char *fileName = NULL;
   /* asprintf(&fileName, "%s-%06d", zLabel, (int)floor(1000*time + 0.5));*/
     asprintf(&fileName, "%s-%06d", zLabel, (int)floor(1000*time ));
    mwplot_plot_gridmap(wp, fileName, nx, ny, A, title);
    if (espera) { wait_for_user(); }
    free(title);
    free(fileName);
  }

void get_field_range(int nx, int ny, double **F, double *vLo, double *vHi, double *vMag)
  {
    (*vLo) = +MAXDOUBLE;
    (*vHi) = -MAXDOUBLE;
    int i, j;
    for (i = 0; i <= nx; i++)
      for (j = 0; j <= ny; j++)
        { double Fij = F[i][j];
          if (! isnan(Fij))
            { if (Fij < (*vLo)) { (*vLo) = Fij; }
              if (Fij > (*vHi)) { (*vHi) = Fij; }
            }
        }
        
    double aLo = fabs(*vLo), aHi = fabs(*vHi);
    (*vMag) = (aLo > aHi ? aLo : aHi);
  }

void get_grid_structure(int nx, int ny, double **F, double **A)
  { /* Copy {F} to {A}, changing NaNs to zeros: */
    int i, j;
    for (i = 0; i <= nx; i++)
      for (j = 0; j <= ny; j++)
        { A[i][j] = (isnan(F[i][j]) ? 0.0 : 1.0); }

    /* DEBUG -- add some {-1} marks to debug plot ranges: */

    auto void mark(int im, int jm);

    void mark(int im, int jm)
      { if (A[im][jm] == 0) { A[im][jm] = -1; } }

    for (i = 0; i < 5; i++)
      { mark(i,0); mark(i,i); mark(0,i);
        mark(nx-i,ny); mark(nx-i,ny-i); mark(nx,ny-i);
      }
  }

void print_grid_statistics(FILE *npf, double t, int nx, int ny, double **A)
  {
    int i, j;
    int nvalid = 0;
    for (i = 0; i <= nx; i++)
      for (j = 0; j <= ny; j++)
        { if (!isnan(A[i][j])) { nvalid = nvalid+1; } }
    fprintf(stderr,  "sparse grid: %7d points\n", nvalid);
    int ntotal=(nx+1)*(ny+1);
    fprintf(stderr,  "full grid:   %7d points\n", ntotal);
    double ratio = ((double)nvalid)/((double)ntotal);
    fprintf(stderr,  "fill ratio:  %7.5f\n", ratio);
    fprintf(stderr, "\n");
    /* fprintf(npf,"%e \t %e \n  %7d points\n", t, ratio, nvalid); */
	for (i = 0; i <= nx; i++)
      for (j = 0; j <= ny; j++)
      fprintf(npf,"%d \t %d \t  %g \n", i, j, A[i][j]);
  }
  
    
void print_dadosEx(int nx, int ny, int k, double **A)
  {
char *name = NULL;
asprintf(&name, "Ex_%04d_%04dc.dat",nx,k);
FILE *p = fopen(name, "w");

//fprintf(p,"# time %g \t  \n",  t);   
int i,j;
double xi, yj;
for (i = 0; i <= nx; i++)
  {xi=(double)i/(double)nx;
  for (j = 0; j <= ny; j++)
    {yj=(double)j/(double)ny;
     fprintf(p,"%g \t %g \t  %g \n", xi, yj, A[i][j]);
    }
  }
fclose(p);      
  }

void print_dadosEy(int nx, int ny, int k, double **A)
  {
char *name = NULL;
asprintf(&name, "Ey_%04d_%04dc.dat",nx,k);
FILE *p = fopen(name, "w");

//fprintf(p,"# time %g \t  \n",  t);   
int i,j;
double xi, yj;
for (i = 0; i <= nx; i++)
  {xi=(double)i/(double)nx;
  for (j = 0; j <= ny; j++)
    {yj=(double)j/(double)ny;
     fprintf(p,"%g \t %g \t  %g \n", xi, yj, A[i][j]);
    }
  }
fclose(p);      
  }

void print_dadosHz(int nx, int ny, int k, double **A)
  {
//p=fopen("dadosHz.dat","w");
char *name = NULL;
asprintf(&name, "Hz_%04d_%04dc.dat",nx,k);
FILE *p = fopen(name, "w");

//fprintf(p,"# time %g \t  \n",  t);   
int i,j;
double xi, yj;
for (i = 0; i <= nx; i++)
  {xi=(double)i/(double)nx;
  for (j = 0; j <= ny; j++)
    {yj=(double)j/(double)ny;
     fprintf(p,"%g \t %g \t  %g \n", xi, yj, A[i][j]);
    }
  }
fclose(p);      
  }