source.c

program FDTD for the photonic crystal structure

#include "./pFDTD.h"

void Gaussian_dipole_source_plane(char *component, int i,int j,float z,float frequency,float phase,long to,long tdecay);
void Lorentzian_dipole_source_plane(char *component, int i,int j,float z,float frequency,float phase,long to,long tdecay);
void Gaussian_dipole_source_line(char *component, float x,int j,float z,float frequency,float phase,long to,long tdecay);
void Lorentzian_dipole_source_line(char *component, float x,int j,float z,float frequency,float phase,long to,long tdecay);

void random_Gaussian_dipole(char *component, float frequency, float tdecay, float x_min, float x_max, float y_min, float y_max, float z_min, float z_max, int gen_number, int seed)
{
	int i;
	float xr, yr, zr, phase;	

	srand(seed); //initialization of random number generator
	
	for(i=0; i<gen_number; i++)
	{
		xr = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(x_max-x_min) + x_min;	
		yr = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(y_max-y_min) + y_min;
		zr = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(z_max-z_min) + z_min;
		phase = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*2*pi;
		Gaussian_dipole_source(component, xr, yr, zr, frequency, phase, 3*tdecay, tdecay);
	}
}

void Gaussian_planewave(char *Ecomp, char *Hcomp, float position, float frequency, long to,long tdecay)
{
	int i, j;
	float phase;	
	
	for(i=0; i<=pisize; i++)
		for(j=0; j<=pjsize; j++)
		{
			phase = 2*pi*(wave_vector_x*i/(lattice_x) + wave_vector_y*j/(lattice_y)); 
			Gaussian_dipole_source_plane(Ecomp, i, j, position, frequency, phase, 3*tdecay, tdecay);
			Gaussian_dipole_source_plane(Hcomp, i, j, position, frequency, phase, 3*tdecay, tdecay);
		}
}

void Gaussian_line_source(char *component, float position_x, float position_z, float frequency, float phase, long to,long tdecay)
{
	int j;
	
	for(j=0; j<=pjsize; j++)
		Gaussian_dipole_source_line(component, position_x, j, position_z, frequency, phase, 3*tdecay, tdecay);
}

void Lorentzian_planewave(char *Ecomp, char *Hcomp, float position, float frequency, long to,long tdecay)
{
	int i, j;
	float phase;	
	
	for(i=0; i<=pisize; i++)
		for(j=0; j<=pjsize; j++)
		{
			phase = 2*pi*(wave_vector_x*i/(lattice_x) + wave_vector_y*j/(lattice_y)); 
			Lorentzian_dipole_source_plane(Ecomp, i, j, position, frequency, phase, to, tdecay);
			Lorentzian_dipole_source_plane(Hcomp, i, j, position, frequency, phase, to, tdecay);
		}
}

void Lorentzian_line_source(char *component, float position_x, float position_z, float frequency, float phase, long to,long tdecay)
{
	int j;
	
	for(j=0; j<=pjsize; j++)
		Lorentzian_dipole_source_line(component, position_x, j, position_z, frequency, phase, to, tdecay);
}

void Gaussian_dipole_source(char *component,float x,float y,float z,float frequency,float phase,long to,long tdecay)
{
	int i,j,k;          // at 3*tdecay Gaussian=0.001 // 

	i=floor(0.5+((x+xcenter)*lattice_x));
	j=floor(0.5+((y+ycenter)*lattice_y));
	k=non_uniform_z_to_i(z);

	if(strcmp(component,"Ex")==0)
	{
		Ex[i][j][k]=Ex[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		Ex[i-1][j][k]=Ex[i-1][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEx[i][j][k]=iEx[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
			iEx[i-1][j][k]=iEx[i-1][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ex")==0)
	{
		Ex[i][j][k]=Ex[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		Ex[i-1][j][k]=Ex[i-1][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEx[i][j][k]=iEx[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
			iEx[i-1][j][k]=iEx[i-1][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		Ey[i][j-1][k]=Ey[i][j-1][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEy[i][j][k]=iEy[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
			iEy[i][j-1][k]=iEy[i][j-1][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		Ey[i][j-1][k]=Ey[i][j-1][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEy[i][j][k]=iEy[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
			iEy[i][j-1][k]=iEy[i][j-1][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		Ez[i][j][k-1]=Ez[i][j][k-1]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEz[i][j][k]=iEz[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
			iEz[i][j][k-1]=iEz[i][j][k-1]+iGauss_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		Ez[i][j][k-1]=Ez[i][j][k-1]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEz[i][j][k]=iEz[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
			iEz[i][j][k-1]=iEz[i][j][k-1]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Hx")==0)
	{
		Hx[i][j][k]=Hx[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hx[i][j-1][k]=Hx[i][j-1][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hx[i][j][k-1]=Hx[i][j][k-1]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hx[i][j-1][k-1]=Hx[i][j-1][k-1]+Gauss_amp(frequency,phase,to,tdecay)/4;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iHx[i][j][k]=iHx[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHx[i][j-1][k]=iHx[i][j-1][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHx[i][j][k-1]=iHx[i][j][k-1]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHx[i][j-1][k-1]=iHx[i][j-1][k-1]+iGauss_amp(frequency,phase,to,tdecay)/4;
		}
	}
	if(strcmp(component,"-Hx")==0)
	{
		Hx[i][j][k]=Hx[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hx[i][j-1][k]=Hx[i][j-1][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hx[i][j][k-1]=Hx[i][j][k-1]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hx[i][j-1][k-1]=Hx[i][j-1][k-1]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iHx[i][j][k]=iHx[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHx[i][j-1][k]=iHx[i][j-1][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHx[i][j][k-1]=iHx[i][j][k-1]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHx[i][j-1][k-1]=iHx[i][j-1][k-1]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
		}
	}
	if(strcmp(component,"Hy")==0)
	{
		Hy[i][j][k]=Hy[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hy[i][j][k-1]=Hy[i][j][k-1]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hy[i-1][j][k]=Hy[i-1][j][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hy[i-1][j][k-1]=Hy[i-1][j][k-1]+Gauss_amp(frequency,phase,to,tdecay)/4;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iHy[i][j][k]=iHy[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHy[i][j][k-1]=iHy[i][j][k-1]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHy[i-1][j][k]=iHy[i-1][j][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHy[i-1][j][k-1]=iHy[i-1][j][k-1]+iGauss_amp(frequency,phase,to,tdecay)/4;
		}
	}
	if(strcmp(component,"-Hy")==0)
	{
		Hy[i][j][k]=Hy[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hy[i][j][k-1]=Hy[i][j][k-1]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hy[i-1][j][k]=Hy[i-1][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hy[i-1][j][k-1]=Hy[i-1][j][k-1]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iHy[i][j][k]=iHy[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHy[i][j][k-1]=iHy[i][j][k-1]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHy[i-1][j][k]=iHy[i-1][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHy[i-1][j][k-1]=iHy[i-1][j][k-1]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
		}
	}
	if(strcmp(component,"Hz")==0)
	{
		Hz[i][j][k]=Hz[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hz[i-1][j][k]=Hz[i-1][j][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hz[i][j-1][k]=Hz[i][j-1][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		Hz[i-1][j-1][k]=Hz[i-1][j-1][k]+Gauss_amp(frequency,phase,to,tdecay)/4;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iHz[i][j][k]=iHz[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHz[i-1][j][k]=iHz[i-1][j][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHz[i][j-1][k]=iHz[i][j-1][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
			iHz[i-1][j-1][k]=iHz[i-1][j-1][k]+iGauss_amp(frequency,phase,to,tdecay)/4;
		}
	}
	if(strcmp(component,"-Hz")==0)
	{
		Hz[i][j][k]=Hz[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hz[i-1][j][k]=Hz[i-1][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hz[i][j-1][k]=Hz[i][j-1][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		Hz[i-1][j-1][k]=Hz[i-1][j-1][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/4;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iHz[i][j][k]=iHz[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHz[i-1][j][k]=iHz[i-1][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHz[i][j-1][k]=iHz[i][j-1][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
			iHz[i-1][j-1][k]=iHz[i-1][j-1][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/4;
		}
	}
}

void Gaussian_dipole_source_plane(char *component, int i,int j,float z,float frequency,float phase,long to,long tdecay)
{
	int k;          // at 3*tdecay Gaussian=0.001 // 

	k=non_uniform_z_to_i(z);

	if(strcmp(component,"Ex")==0)
	{
		Ex[i][j][k]=Ex[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if(i>=1)
			Ex[i-1][j][k]=Ex[i-1][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEx[i][j][k]=iEx[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
			if(i>=1)
				iEx[i-1][j][k]=iEx[i-1][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ex")==0)
	{
		Ex[i][j][k]=Ex[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if(i>=1)
			Ex[i-1][j][k]=Ex[i-1][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEx[i][j][k]=iEx[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
			if(i>=1)
				iEx[i-1][j][k]=iEx[i-1][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if(j>=1)
			Ey[i][j-1][k]=Ey[i][j-1][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEy[i][j][k]=iEy[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
			if(j>=1)
				iEy[i][j-1][k]=iEy[i][j-1][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if(j>=1)
			Ey[i][j-1][k]=Ey[i][j-1][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEy[i][j][k]=iEy[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
			if(j>=1)
				iEy[i][j-1][k]=iEy[i][j-1][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if(k>=1)
			Ez[i][j][k-1]=Ez[i][j][k-1]+Gauss_amp(frequency,phase,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEz[i][j][k]=iEz[i][j][k]+iGauss_amp(frequency,phase,to,tdecay)/2;
			if(k>=1)
				iEz[i][j][k-1]=iEz[i][j][k-1]+iGauss_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if(k>=1)
			Ez[i][j][k-1]=Ez[i][j][k-1]+Gauss_amp(frequency,phase+pi,to,tdecay)/2;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{
			iEz[i][j][k]=iEz[i][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
			if(k>=1)
				iEz[i][j][k-1]=iEz[i][j][k-1]+iGauss_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Hx")==0)
	{