source.c

program FDTD for the photonic crystal structure

		if(i>=1)
			Hz[i-1][j][k]=Hz[i-1][j][k]+Gauss_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
		if(j>=1)
			Hz[i][j-1][k]=Hz[i][j-1][k]+Gauss_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
		if(i>=1 && j>=1)
			Hz[i-1][j-1][k]=Hz[i-1][j-1][k]+Gauss_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/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)*sqrt(eo/uo)/4;
			if(i>=1)
				iHz[i-1][j][k]=iHz[i-1][j][k]+iGauss_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
			if(j>=1)
				iHz[i][j-1][k]=iHz[i][j-1][k]+iGauss_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
			if(i>=1 && j>=1)
				iHz[i-1][j-1][k]=iHz[i-1][j-1][k]+iGauss_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
		}
	}
}

void Lorentzian_dipole_source(char *component,float x,float y,float z,float frequency,float phase,long to,long tdecay)
{
	int i,j,k;         //at 8*tdecay Lorentzian=0.0001 // 

	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]+Lorentz_amp(frequency,phase,to,tdecay)/2;
		Ex[i-1][j][k]=Ex[i-1][j][k]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/2;
			iEx[i-1][j][k]=iEx[i-1][j][k]+iLorentz_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ex")==0)
	{
		Ex[i][j][k]=Ex[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/2;
		Ex[i-1][j][k]=Ex[i-1][j][k]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
			iEx[i-1][j][k]=iEx[i-1][j][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/2;
		Ey[i][j-1][k]=Ey[i][j-1][k]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/2;
			iEy[i][j-1][k]=iEy[i][j-1][k]+iLorentz_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/2;
		Ey[i][j-1][k]=Ey[i][j-1][k]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
			iEy[i][j-1][k]=iEy[i][j-1][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/2;
		Ez[i][j][k-1]=Ez[i][j][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/2;
			iEz[i][j][k-1]=iEz[i][j][k-1]+iLorentz_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/2;
		Ez[i][j][k-1]=Ez[i][j][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
			iEz[i][j][k-1]=iEz[i][j][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Hx")==0)
	{
		Hx[i][j][k]=Hx[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hx[i][j-1][k]=Hx[i][j-1][k]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hx[i][j][k-1]=Hx[i][j][k-1]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hx[i][j-1][k-1]=Hx[i][j-1][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHx[i][j-1][k]=iHx[i][j-1][k]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHx[i][j][k-1]=iHx[i][j][k-1]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHx[i][j-1][k-1]=iHx[i][j-1][k-1]+iLorentz_amp(frequency,phase,to,tdecay)/4;
		}
	}
	if(strcmp(component,"-Hx")==0)
	{
		Hx[i][j][k]=Hx[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hx[i][j-1][k]=Hx[i][j-1][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hx[i][j][k-1]=Hx[i][j][k-1]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hx[i][j-1][k-1]=Hx[i][j-1][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHx[i][j-1][k]=iHx[i][j-1][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHx[i][j][k-1]=iHx[i][j][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHx[i][j-1][k-1]=iHx[i][j-1][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
		}
	}
	if(strcmp(component,"Hy")==0)
	{
		Hy[i][j][k]=Hy[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hy[i][j][k-1]=Hy[i][j][k-1]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hy[i-1][j][k]=Hy[i-1][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hy[i-1][j][k-1]=Hy[i-1][j][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHy[i][j][k-1]=iHy[i][j][k-1]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHy[i-1][j][k]=iHy[i-1][j][k]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHy[i-1][j][k-1]=iHy[i-1][j][k-1]+iLorentz_amp(frequency,phase,to,tdecay)/4;
		}
	}
	if(strcmp(component,"-Hy")==0)
	{
		Hy[i][j][k]=Hy[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hy[i][j][k-1]=Hy[i][j][k-1]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hy[i-1][j][k]=Hy[i-1][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hy[i-1][j][k-1]=Hy[i-1][j][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHy[i][j][k-1]=iHy[i][j][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHy[i-1][j][k]=iHy[i-1][j][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHy[i-1][j][k-1]=iHy[i-1][j][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
		}
	}
	if(strcmp(component,"Hz")==0)
	{
		Hz[i][j][k]=Hz[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hz[i-1][j][k]=Hz[i-1][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hz[i][j-1][k]=Hz[i][j-1][k]+Lorentz_amp(frequency,phase,to,tdecay)/4;
		Hz[i-1][j-1][k]=Hz[i-1][j-1][k]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHz[i-1][j][k]=iHz[i-1][j][k]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHz[i][j-1][k]=iHz[i][j-1][k]+iLorentz_amp(frequency,phase,to,tdecay)/4;
			iHz[i-1][j-1][k]=iHz[i-1][j-1][k]+iLorentz_amp(frequency,phase,to,tdecay)/4;
		}
	}
	if(strcmp(component,"-Hz")==0)
	{
		Hz[i][j][k]=Hz[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hz[i-1][j][k]=Hz[i-1][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hz[i][j-1][k]=Hz[i][j-1][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/4;
		Hz[i-1][j-1][k]=Hz[i-1][j-1][k]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHz[i-1][j][k]=iHz[i-1][j][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHz[i][j-1][k]=iHz[i][j-1][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
			iHz[i-1][j-1][k]=iHz[i-1][j-1][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/4;
		}
	}
}

void Lorentzian_dipole_source_plane(char *component,int i,int j,float z,float frequency,float phase,long to,long tdecay)
{
	int k;         //at 8*tdecay Lorentzian=0.0001 // 

	k=non_uniform_z_to_i(z);

	if(strcmp(component,"Ex")==0)
	{
		Ex[i][j][k]=Ex[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/2;
		if(i>=1)
			Ex[i-1][j][k]=Ex[i-1][j][k]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/2;
			if(i>=1)
				iEx[i-1][j][k]=iEx[i-1][j][k]+iLorentz_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ex")==0)
	{
		Ex[i][j][k]=Ex[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/2;
		if(i>=1)
			Ex[i-1][j][k]=Ex[i-1][j][k]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
			if(i>=1)
				iEx[i-1][j][k]=iEx[i-1][j][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/2;
		if(j>=1)
			Ey[i][j-1][k]=Ey[i][j-1][k]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/2;
			if(j>=1)
				iEy[i][j-1][k]=iEy[i][j-1][k]+iLorentz_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ey")==0)
	{
		Ey[i][j][k]=Ey[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/2;
		if(j>=1)
			Ey[i][j-1][k]=Ey[i][j-1][k]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
			if(j>=1)
				iEy[i][j-1][k]=iEy[i][j-1][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)/2;
		if(k>=1)
			Ez[i][j][k-1]=Ez[i][j][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase,to,tdecay)/2;
			if(k>=1)
				iEz[i][j][k-1]=iEz[i][j][k-1]+iLorentz_amp(frequency,phase,to,tdecay)/2;
		}
	}
	if(strcmp(component,"-Ez")==0)
	{
		Ez[i][j][k]=Ez[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)/2;
		if(k>=1)
			Ez[i][j][k-1]=Ez[i][j][k-1]+Lorentz_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]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
			if(k>=1)
				iEz[i][j][k-1]=iEz[i][j][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)/2;
		}
	}
	if(strcmp(component,"Hx")==0)
	{
		Hx[i][j][k]=Hx[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		if(j>=1)
			Hx[i][j-1][k]=Hx[i][j-1][k]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		if(k>=1)
			Hx[i][j][k-1]=Hx[i][j][k-1]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		if(j>=1 && k>=1)
			Hx[i][j-1][k-1]=Hx[i][j-1][k-1]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/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]+iLorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
			if(j>=1)
				iHx[i][j-1][k]=iHx[i][j-1][k]+iLorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
			if(k>=1)
				iHx[i][j][k-1]=iHx[i][j][k-1]+iLorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
			if(j>=1 && k>=1)
				iHx[i][j-1][k-1]=iHx[i][j-1][k-1]+iLorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		}
	}
	if(strcmp(component,"-Hx")==0)
	{
		Hx[i][j][k]=Hx[i][j][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
		if(j>=1)
			Hx[i][j-1][k]=Hx[i][j-1][k]+Lorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
		if(k>=1)
			Hx[i][j][k-1]=Hx[i][j][k-1]+Lorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
		if(j>=1 && k>=1)
			Hx[i][j-1][k-1]=Hx[i][j-1][k-1]+Lorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/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]+iLorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
			if(j>=1)
				iHx[i][j-1][k]=iHx[i][j-1][k]+iLorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
			if(k>=1)
				iHx[i][j][k-1]=iHx[i][j][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
			if(j>=1 && k>=1)
				iHx[i][j-1][k-1]=iHx[i][j-1][k-1]+iLorentz_amp(frequency,phase+pi,to,tdecay)*sqrt(eo/uo)/4;
		}
	}
	if(strcmp(component,"Hy")==0)
	{
		Hy[i][j][k]=Hy[i][j][k]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		if(k>=1)
			Hy[i][j][k-1]=Hy[i][j][k-1]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		if(i>=1)
			Hy[i-1][j][k]=Hy[i-1][j][k]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		if(i>=1 && k>=1)
			Hy[i-1][j][k-1]=Hy[i-1][j][k-1]+Lorentz_amp(frequency,phase,to,tdecay)*sqrt(eo/uo)/4;
		if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		{