/* Last edited on 2006-06-08 by anamaria */ 


/*
 * Integrates the UMPL system of equations  in 2D
 * with Gaussian initial condition 
 * 14/05/2006
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "GradFunc.h"

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

int main()
{
  /* Fundamental constants */
  double cc = 2.99792458e8;                 /*speed of light in free space */
  double lnR0 = log(1.0e-7);                /* ln(R(0)) */
  double mu0 = 4.0 * M_PI * 1.0e-7;         /* permeability of free space */
  double eps0 = 1.0 / (cc*cc*mu0);          /* permittivity of free  space */ 
  double eta = sqrt(mu0/eps0);              /* wave impedance in free space ???*/

  double mur  = 1.0;                        /* relative permeability in the PML region*/
  double epsr = 1.0 ;                       /* relative permittivity in the PML region*/

  int n = 64;                               /* the matriz size nxn*/
  int l = 10;                                /* PML size = l*dx */  
  double d = (double) l/(double) n;
  double or = 4.00;			    /* order of PML; Taflove is m*/

  int interp_points = 4;                    /* 4 point interpolation ( cubic) */
  double dx = 1.0/(n);				/* space increment in x-direction */
  double dy = 1.0/(n);				/* space increment in y-direction */
  double dt = 0.707* dx / cc;			/* time step */
  int tmax =(int)(8.0e-9 / dt);			/* total number of time steps */

 
  double kappamax = 1.0;                             /* maximum value for kappa at the PML boundary*/
  double sigmax = (or + 1.0)/(150.0*M_PI*sqrt(epsr)*dx);  /* maximum value for sigma at the PML boundary*/

 
  int i, j, k;	
  
 		
  double **DDx=matrix_alloc(n+1, n+1);              /* memoria para a matriz do Dx*/
  double **DDy=matrix_alloc(n+1, n+1);              /* idem para o Dy */
  double **BBz=matrix_alloc(n+1, n+1);              /* idem para o Bz*/
  
  grade_t *Hz = cria_grade(n, l, mu0,  BBz);   /* Componente Z do campo H. */
  grade_t *Ex = cria_grade(n, l, eps0, DDx);   /* Componente X do campo E. */
  grade_t *Ey = cria_grade(n, l, eps0, DDy);   /* Componente Y do campo E. */

  double **temp=matrix_alloc(n+1, n+1);
  double **tplt=matrix_alloc(n+1, n+1);

  /*definindo as variaveis da regiao de PML*/
  double sig[l];
  double kappa[l];
  double ck1[n+1];
  double ck2[n+1];

  
  double dyhz;
  double dxhz;
  double dyex;
  double dxey;

  /* Set up the PML material coefficients */
  for (i=0; i <l; i++)
    {
      sig[i] = sigmax*pow(((l-i)*1.0)/(l),or);
      kappa[i] = 1.0+(kappamax-1.0)*pow(((l-i)*1.0)/(l),or);
    }

  for (i=0; i <l; i++)
    { ck1[i]=kappa[i]-sig[i]*dt/(2*eps0);
      ck1[n-i]=ck1[i];
      ck2[i]=kappa[i]+sig[i]*dt/(2*eps0);
      ck2[n-i]=ck2[i];
    }

  for (i=l; i <=n-l; i++) ck1[i]=ck2[i]=1.0;

  /* Inicializando as variáveis BBz, DDx e DDy*/
  for (i=0; i <=n; i++)
    for (j=0; j <=n; j++)
      {
        BBz[i][j]= mu0*exp(-(i*dx-0.5)*(i*dx-0.5)*150)*exp(-(j*dy-0.5)*(j*dy-0.5)*150);
        DDx[i][j]=0.0;
        DDy[i][j]=0.0;
      }

  /* Inicializando as variáveis na regiao de UPML*/
  for (i=0; i <= n; i++)
    for (j=0; j <=n; j++)
      { if ((i < l) || (i > n-l) || (j < l) || (j > n-l))
          { /* [i][j] está na borda UPML */
            mset(Hz, i, j, BBz[i][j]/mu0);
            mset(Ex, i, j, DDx[i][j]/eps0);
            mset(Ey, i, j, DDy[i][j]/eps0);
          }
      }

  /*Graficando a condição inicial do campo magnetico*/
  FILE *pp = abre_grafico(n, dx, dy);
    mostra_grade(pp, "Hz", 0.0, n, dx, dy, 1.0, Hz, tplt, 1);
    	
  /* BEGIN TIME-STEPPING LOOP */
  int debug = 1;
  
  for (k=1; k<= tmax; k++)
    {
      double tempo = k * dt * 1e9;
      /* int esp = (tempo > 1.2); */
      int esp = (tempo > 3.0);

      /*fazendo o update do BBz*/
      for (i=0; i <=n; i++)
        for (j=0; j <=n; j++)
          { temp[i][j]= BBz[i][j];
            dyex= DerivAntSimet(Ex, dy, i, j, 0); 
            dxey= DerivAntSimet(Ey, dx, i, j, 1);
            BBz[i][j]=ck1[i]*temp[i][j]/ck2[i]+dt*(dyex-dxey)/ck2[i];
          }
      if (debug) mostra_matriz(pp, "BBz", tempo, n, dx, dy, 2.0e-6, BBz, esp);

      /*fazendo o update do Hz */
      for (i=0; i <=n; i++)
        for (j=0; j <=n; j++)
          { if ((i < l) || (i > n-l) || (j < l) || (j > n-l))
              { /* [i][j] está na borda. */
                double Hz1 = mget(Hz, i, j);
                double Hz2 = ck1[j]*Hz1/ck2[j] + (BBz[i][j]-temp[i][j])/(ck2[j]*mu0*mur);
                mset(Hz, i, j, Hz2);
              }
          }
     if (debug) mostra_grade(pp, "Hz", tempo, n, dx, dy, 1.0, Hz, tplt, esp);
       
      /*fazendo o update do DDx*/
      for (i=0; i <=n; i++)
        for (j=0; j <=n; j++)
          {temp[i][j]=DDx[i][j];
            dyhz=DerivSimet(Hz,  dy, i, j, 0);
            DDx[i][j]=ck1[j]*temp[i][j]/ck2[j] + dt*dyhz/ck2[j];
	    
          }
      if (debug) mostra_matriz(pp, "DDx", tempo, n, dx, dy, 2.0e-9, DDx, esp);

      /*fazendo o update do Ex */
      for (i=0; i <=n; i++)
        for (j=0; j <=n; j++)
          { if ((i < l) || (i > n-l) || (j < l) || (j > n-l))
              { /* [i][j] está na borda. */
                double Ex1 = mget(Ex, i, j);
                double Ex2 = Ex1 + (ck2[i]*DDx[i][j]- ck1[i]*temp[i][j])/(eps0*epsr);
                mset(Ex, i, j, Ex2);
              }
          }
      if (debug) mostra_grade(pp, "Ex", tempo, n, dx, dy, 100.0, Ex, tplt, esp);

      /*fazendo o update do DDy*/
	 
      for (i=0; i <=n; i++)
        for (j=0; j <=n; j++)
          { temp[i][j]=DDy[i][j];
            dxhz=DerivSimet(Hz,  dx, i, j, 1);
            DDy[i][j]=temp[i][j] - dt*dxhz;
          }
      if (debug) mostra_matriz(pp, "DDy", tempo, n, dx, dy, 2.0e-9, DDy, esp);

      /*fazendo o update do Ey */
      for (i=0; i <=n; i++)
        for (j=0; j <=n; j++)
          { if ((i < l) || (i > n-l) || (j < l) || (j > n-l))
              { /* [i][j] está na borda. */
                double Ey1 = mget(Ey, i, j);
                double Ey2 = ck1[i]*Ey1/ck2[i]+(ck2[j]*DDy[i][j]- ck1[j]*temp[i][j])/(ck2[i]*eps0*epsr);
                mset(Ey, i, j, Ey2);
              }
          }
      if (debug) mostra_grade(pp, "Ey", tempo, n, dx, dy, 100.0, Ey, tplt, esp);

      mostra_grade(pp, "Hz", tempo, n, dx, dy, 1.0, Hz, tplt, esp);

    } /* fim da evolucao no tempo */	
	
  pclose(pp);

  return 0;
}