// The following code was created by Rolf Timp in June 2007 for AP273 
// at Stanford as a demonstration of the utility of finite difference 
// time domain methods in calculating 2D photonic band gaps. The code is 
// written in C, and was compiled using Dev-C++ - a C and C++ compiler. 
// Dev-C++ is free software distributed under the GNU General Public 
// License. It can be obtained at 
// http://www.bloodshed.net/dev/devcpp.html

// Note: Compiles under g++ but not gcc. - RBL 25 Jun 07.

// The code uses three different arrays to store three different fields: 
// Ez, Hx, and Hy. Using a discretized form of Maxwell's equations, an 
// electric field delta function is introduced into the system, and the 
// system's impulse response is measured. Currently, the program creates 
// a file called Efield.txt which stores the final state of the electric 
// field after the specified number of timesteps. However, one could 
// also save the state of the electric field at every time step which 
// would permit one to actually view the delta function's propagation.

#include
#include
#include

main()
{

// Constants
float pi = 3.14159;

float c = 3*100000000;
float mu = 4*pi*0.0000001;
float e0 = 8.85*0.000000000001; // correct exponent = 0.00000000001

// Numbers of Timesteps
float tlength = 100;

// Simulation's Resolution
float deltax2 = 2*0.00000001;
float deltay2 = 2*0.00000001;
float deltat2 = deltax2/(1.41421356*c);

// Size of the Crystal in units of deltax2
float xlength2 = 201;
float ylength2 = 201;

// Grids storing the magnetic and electric fields on the gridpoints
float Hx2[201][201] = {0}; // -1 on y side
float Hy2[201][201] = {0}; // -1 on x side
float Ez2[201][201] = {0};
float Ezx2[201][201] = {0};
float Ezy2[201][201] = {0};

float epsilon[201][201] = {e0}; // the dielectric constant at each point on the grid	
float sigmax[201][201] = {0}; // a constant used for implementation of the perfectly conducting layer boundary condition (PML)
float sigmay[201][201] = {0}; // used for PML
float sigmax2[201][201] = {0}; // used for PML
float sigmay2[201][201] = {0}; // used for PML

float Ezxcoeff1, Ezxcoeff2, Ezycoeff1, Ezycoeff2, Hycoeff1, Hycoeff2, Hxcoeff1, Hxcoeff2;

int i, j, n;

// Creating the periodic variations in the photonic crystal's dielectric constant
for ( j = 0 ; j <= ylength2-1; j++){
    for(i = 0; i <= xlength2-1; i++){  
        
        if(((i%29 - 10)*(i%29 - 10) + (j%29 - 10)*(j%29 - 10)) < 33){
            epsilon[i][j] = 8.9*e0;}
        else
            epsilon[i][j] = e0;
            
            }
        }    

// The simulation itself, based off a discretization of Maxwell's equations.

printf("The simulation is running...\n");

for (n = 1; n<= tlength-1; n++){
printf("Time step #:%d\n", n);

    // The magnetic field must be updated first, using the discretization I chose.
    for (j = 0; j<=ylength2-1; j++){
    	for (i = 0; i<=xlength2-1; i++){ 

        int  j_retarded = j-1;
        int  i_retarded = i-1;

        // coefficients used to put Maxwell's equations in a neater, more recognizable form

        Hycoeff1 = (2*mu-deltat2*sigmax2[i][j])/(2*mu+deltat2*sigmax2[i][j]);
        Hycoeff2 = (2*deltat2)/(deltax2*(2*mu+deltat2*sigmax2[i][j]));
        Hxcoeff1 = (2*mu-deltat2*sigmay2[i][j])/(2*mu+deltat2*sigmay2[i][j]);
        Hxcoeff2 = (2*deltat2)/(deltay2*(2*mu+deltat2*sigmay2[i][j]));
        
        // The update equation for the x-component of the magnetic field
        if(j>=1&&i<=xlength2-2)
          Hx2[i][j] = Hxcoeff1*Hx2[i][j] - Hxcoeff2*(Ez2[i][j+1] - Ez2[i][j]); 
    
        // The update equation for the y-component of the magnetic field
        if(i>=1&&j<=ylength2-2)
            Hy2[i][j] = Hycoeff1*Hy2[i][j] + Hycoeff2*(Ez2[i+1][j] - Ez2[i][j]);

      }
    }
    
    for (j = 0; j<=ylength2-1; j++){
        for (i = 0; i<=xlength2-1; i++){ 
            
                    int  j_retarded = j-1;
                    int  i_retarded = i-1;

        // coefficients used to put Maxwell's equations in a neater, more recognizable form            
            Ezxcoeff1 = ((2*epsilon[i][j] - deltat2*sigmax[i][j])/(2*epsilon[i][j] + deltat2*sigmax[i][j]));
            Ezxcoeff2 = ((2*deltat2)/(deltax2*(2*epsilon[i][j] + deltat2*sigmax[i][j])));
            Ezycoeff1 = ((2*epsilon[i][j] - deltat2*sigmay[i][j])/(2*epsilon[i][j] + deltat2*sigmay[i][j]));;
            Ezycoeff2 = ((2*deltat2)/(deltay2*(2*epsilon[i][j] + deltat2*sigmay[i][j])));
  
        //  E-field update equations          
            if(i>=1 && i < 201 && j>=1 && j < 201){
                Ezx2[i][j] = Ezxcoeff1*Ezx2[i][j] + Ezxcoeff2*(Hy2[i][j] - Hy2[i_retarded][j]);
                Ezy2[i][j] = Ezycoeff1*Ezy2[i][j] - Ezycoeff2*((Hx2[i][j] - Hx2[i][j_retarded]));
                Ez2[i][j]  =  Ezx2[i][j] + Ezy2[i][j]; 
            }
  
        }
    }
    
   	// The spatial and temporal delta function which is used to determine the system's frequency response.
   	// It's the system's input.
   	
    float x = (0.2*(n-30));
    Ez2[107][94] = 100*exp(-x*x); //100*(1+x+x*x/2+x*x*x/6+x*x*x*x/24+x*x*x*x*x/120+x*x*x*x*x*x/720+x*x*x*x*x*x*x/(7*720)+x*x*x*x*x*x*x*x/(8*7*720)+x*x*x*x*x*x*x*x*x/(9*8*7*720)+x*x*x*x*x*x*x*x*x*x/(10*9*8*7*720));  
    
}

// Writing the Efield to the file Efield.txt

    FILE *file;
    file = fopen("Efield.txt","a+");
    
        for (j = 0; j<=ylength2-1; j++){
        	for (i = 0; i<=xlength2-1; i++){
            fprintf(file,"%e ",Ez2[i][j]," "); 
            }
            fprintf(file,"\n");
         }
         
         fclose(file);
     
}