output.c

program FDTD for the photonic crystal structure

#include "./pFDTD.h"

float efieldx(int i,int j,int k);
float efieldy(int i,int j,int k);
float efieldz(int i,int j,int k);
float jfieldx(int i,int j,int k);
float jfieldy(int i,int j,int k);
float jfieldz(int i,int j,int k);
float hfieldx(int i,int j,int k);
float hfieldy(int i,int j,int k);
float hfieldz(int i,int j,int k);
float iefieldx(int i,int j,int k);
float iefieldy(int i,int j,int k);
float iefieldz(int i,int j,int k);
float ijfieldx(int i,int j,int k);
float ijfieldy(int i,int j,int k);
float ijfieldz(int i,int j,int k);
float ihfieldx(int i,int j,int k);
float ihfieldy(int i,int j,int k);
float ihfieldz(int i,int j,int k);
float eps(int i,int j,int k);
float meps(int i,int j,int k); //For Hz_parity()

float grid_value_periodic_x(char *component,int i,int j,int k, int h);
float grid_value_periodic_y(char *component,int i,int j,int k, int h);
float grid_value_periodic_xy(char *component,int i,int j,int k, int h, int v);

int find_k_projection(int i, int j, int k_range);

void real_space_param(float a_nm, float w_n)
{
	// a_nm is in the unit of 'nanometer'
	// w_n is the normalized frequency
	FILE *PARAM;

	PARAM=fopen("Real_Space_Param.dat","wt");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The FDTD calculation domain \n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"Lx =%lf (a) \n",xsize);
	fprintf(PARAM,"Ly =%lf (a) \n",ysize);
	fprintf(PARAM,"Lz =%lf (a) \n",zsize);
	fprintf(PARAM,"---------------------------------\n");
	fprintf(PARAM,"Lx =%lf (nm) \n",xsize*a_nm);
	fprintf(PARAM,"Ly =%lf (nm) \n",ysize*a_nm);
	fprintf(PARAM,"Lz =%lf (nm) \n",zsize*a_nm);
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The PML size \n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"pmlil =%d (u), pmlir =%d (u) \n",pmlil,pmlir);
	fprintf(PARAM,"pmljl =%d (u), pmljr =%d (u) \n",pmljl,pmljr);
	fprintf(PARAM,"pmlkl =%d (u), pmlkr =%d (u) \n",pmlkl,pmlkr);
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The speed of light \n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"c =%E (m/s) \n",light_speed);
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The Stability parameter 'S' \n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"S =%lf (1/q) \n",S_factor);
	fprintf(PARAM,"S =%E (1/m) \n",S_factor*ds_x*lattice_x/(a_nm*1E-9));
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The size of the FDTD grid \n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"dx =%lf (q) \n",ds_x);
	fprintf(PARAM,"dy =%lf (q) \n",ds_y);
	fprintf(PARAM,"dz =%lf (q) \n",ds_z);
	fprintf(PARAM,"---------------------------------\n");
	fprintf(PARAM,"dx =%lf (nm) \n",a_nm/lattice_x);
	fprintf(PARAM,"dy =%lf (nm) \n",a_nm/lattice_y);
	fprintf(PARAM,"dz =%lf (nm) \n",a_nm/lattice_nz);
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The size of the FDTD time step \n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"dt =%g (sec) \n",(a_nm*1E-9)/(light_speed*S_factor*ds_x*lattice_x));
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The dipole wavelength in vacuum (in FDTD)\n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"lambda =%lf (u) \n",lattice_x/w_n);
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The dipole frequency (in FDTD)\n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"f =%lf (1/update) \n",w_n/(S_factor*ds_x*lattice_x));
	fprintf(PARAM,"T =%lf (update) \n",(S_factor*ds_x*lattice_x)/w_n);
	fprintf(PARAM,"\n");

	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"The Origin FFT correction factor \n");
	fprintf(PARAM,"=================================\n");	
	fprintf(PARAM,"x%lf\n",S_factor*ds_x*lattice_x);
	fprintf(PARAM,"\n");

	fclose(PARAM);
}

int get_period_in_update(float w_n)
{
	return((int)(S_factor*ds_x*lattice_x)/w_n);	
}

void out_epsilon(char *plane,float value,char *name)
{
	FILE *file;
	int i,j,k;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	if(use_periodic_x == 1)
		i_range = isize;
	if(use_periodic_y == 1)
		j_range = jsize; 

	file=fopen(name,"w");

	if(strcmp("x",plane)==0)
	{
		i=floor(0.5+((value+xcenter)*lattice_x));
		for(k=k_range;k>=1;k--)
		{
			for(j=1;j<=j_range;j++)
			{
				if(meps(i,j,k)==0.0)
					fprintf(file,"%g ",eps(i,j,k));
				else
					fprintf(file,"%g ",meps(i,j,k));
			}
			fprintf(file,"\n");
		}
	}
	if(strcmp("y",plane)==0)
	{
		j=floor(0.5+((value+ycenter)*lattice_y));
		for(k=k_range;k>=1;k--)
		{
			for(i=1;i<=i_range;i++)
			{
				if(meps(i,j,k)==0.0)
					fprintf(file,"%g ",eps(i,j,k));
				else
					fprintf(file,"%g ",meps(i,j,k));
			}
			fprintf(file,"\n");
		}
	}

	if(strcmp("z",plane)==0)
	{
		k=non_uniform_z_to_i(value);
		for(j=j_range;j>=1;j--)
		{
			for(i=1;i<=i_range;i++)
			{
				if(meps(i,j,k)==0.0)
					fprintf(file,"%g ",eps(i,j,k));
				else
					fprintf(file,"%g ",meps(i,j,k));
			}
			fprintf(file,"\n");
		}
	}
	fclose(file);
	printf("out_epsilon...ok\n");
}

void out_epsilon_periodic(char *plane,float value,char *name, int m_h, int m_v)
{
	FILE *file;
	int i,j,k;
	int v,h;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	file=fopen(name,"w");

	if(strcmp("x",plane)==0)
	{
		i=floor(0.5+((value+xcenter)*lattice_x));
		for(v=0; v<m_v; v++)
		{
			for(k=k_range;k>=1;k--)
			{
				for(h=0; h<m_h; h++)
				{
					for(j=1;j<=j_range;j++)
					{
						if(meps(i,j,k)==0.0)
							fprintf(file,"%g ",eps(i,j,k));
						else
							fprintf(file,"%g ",meps(i,j,k));
					}
				}
				fprintf(file,"\n");
			}
		}
	}

	if(strcmp("y",plane)==0)
	{
		j=floor(0.5+((value+ycenter)*lattice_y));
		for(v=0; v<m_v; v++)
		{
			for(k=k_range;k>=1;k--)
			{
				for(h=0; h<m_h; h++)
				{
					for(i=1;i<=i_range;i++)
					{
						if(meps(i,j,k)==0.0)
							fprintf(file,"%g ",eps(i,j,k));
						else
							fprintf(file,"%g ",meps(i,j,k));
					}
				}
				fprintf(file,"\n");
			}
		}
	}

	if(strcmp("z",plane)==0)
	{
		k=non_uniform_z_to_i(value);
		for(v=0; v<m_v; v++)
		{
			for(j=j_range;j>=1;j--)
			{
				for(h=0; h<m_h; h++)
				{
					for(i=1;i<=i_range;i++)
					{
						if(meps(i,j,k)==0.0)
							fprintf(file,"%g ",eps(i,j,k));
						else
							fprintf(file,"%g ",meps(i,j,k));
					}
				}
				fprintf(file,"\n");
			}
		}
	}
	fclose(file);
	printf("out_epsilon...ok\n");
}

void out_epsilon_projection(char *dirc, char *name)
{
	FILE *file;
	int i,j,k;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	if(use_periodic_x == 1)
		i_range = isize;
	if(use_periodic_y == 1)
		j_range = jsize; 

	file=fopen(name,"w");

	if(strcmp("+z",dirc)==0)
	{
		for(j=j_range;j>=1;j--)
		{
			for(i=1;i<=i_range;i++)
			{
				k= find_k_projection(i,j,k_range);
				if(meps(i,j,k)==0.0)
					fprintf(file,"%g ",eps(i,j,k));
				else
					fprintf(file,"%g ",meps(i,j,k));
			}
			fprintf(file,"\n");
		}
	}
	fclose(file);
	printf("out_epsilon...ok\n");
}

void out_plane(char *component,char *plane,float value,char *lastname)
{
	FILE *file;
	char name[20];
	int i,j,k;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	if(use_periodic_x == 1)
		i_range = isize;
	if(use_periodic_y == 1)
		j_range = jsize; 
	
	sprintf(name,"%07d",t);
	strcat(name,lastname);
	file=fopen(name,"w");

	if(strcmp("x",plane)==0)
	{
		i=floor(0.5+((value+xcenter)*lattice_x));
		for(k=k_range;k>=1;k--)
		{
			for(j=1;j<=j_range;j++)	fprintf(file,"%g ",grid_value(component,i,j,k));
			fprintf(file,"\n");
		}
	}

	if(strcmp("y",plane)==0)
	{
		j=floor(0.5+((value+ycenter)*lattice_y));
		for(k=k_range;k>=1;k--)
		{
			for(i=1;i<=i_range;i++)	fprintf(file,"%g ",grid_value(component,i,j,k));
			fprintf(file,"\n");
		}
	}

	if(strcmp("z",plane)==0)
	{
		k=non_uniform_z_to_i(value);
		for(j=j_range;j>=1;j--)
		{
			for(i=1;i<=i_range;i++) fprintf(file,"%g ",grid_value(component,i,j,k));	
			fprintf(file,"\n");
		}
	}
	fclose(file);
	printf("out %s...ok\n",component);
}

void out_plane_periodic(char *component,char *plane,float value,char *lastname, int m_h, int m_v)
{
	FILE *file;
	char name[20];
	int i,j,k;
	int h, v;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	sprintf(name,"%07d",t);
	strcat(name,lastname);
	file=fopen(name,"w");

	if(strcmp("x",plane)==0)
	{
		i=floor(0.5+((value+xcenter)*lattice_x));
		for(v=0; v<1; v++)  // no periodic for z-direction
		{
			for(k=k_range;k>=1;k--)
			{
				for(h=0; h<m_h; h++)
				{
					for(j=1;j<=j_range;j++)
						fprintf(file,"%g ",grid_value_periodic_y(component,i,j,k,h));	
				}
				fprintf(file,"\n");
			}
		}
	}

	if(strcmp("y",plane)==0)
	{
		j=floor(0.5+((value+ycenter)*lattice_y));
		for(v=0; v<1; v++)  // no periodic for z-direction
		{
			for(k=k_range;k>=1;k--)
			{
				for(h=0; h<m_h; h++)
				{
					for(i=1;i<=i_range;i++)
						fprintf(file,"%g ",grid_value_periodic_x(component,i,j,k,h));	
				}
				fprintf(file,"\n");
			}
		}
	}

	if(strcmp("z",plane)==0)
	{
		k=non_uniform_z_to_i(value);
		for(v=0; v<m_v; v++)
		{
			for(j=j_range;j>=1;j--)
			{
				for(h=0; h<m_h; h++)
				{
					for(i=1;i<=i_range;i++)
						fprintf(file,"%g ",grid_value_periodic_xy(component,i,j,k,h,v));	
				}
				fprintf(file,"\n");
			}
		}
	}
	fclose(file);
	printf("out %s...ok\n",component);
}

void out_plane_time_average(char *component,char *plane,float value, long start, long end, float ***field_avg, char *lastname)
{
	FILE *file;
	char name[20];
	int i,j,k;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	if(use_periodic_x == 1)
		i_range = isize;
	if(use_periodic_y == 1)
		j_range = jsize; 

	if(strcmp("x",plane)==0)
	{
		if(t == start)
		{
			///// *field_avg[j][k] /////
			(*field_avg) = (float **)malloc(sizeof(float *)*(j_range+1));
			for(j=0; j<j_range+1; j++)
				(*field_avg)[j] = (float *)malloc(sizeof(float)*(k_range+1));

			for(j=0; j<j_range+1; j++) 
				for(k=0; k<k_range+1; k++)
					(*field_avg)[j][k] = 0.0;
		}

		if(t>start && t<=end)
		{
			i=floor(0.5+((value+xcenter)*lattice_x));
			for(k=k_range;k>=1;k--)
			{
				for(j=1;j<=j_range;j++)  
					(*field_avg)[j][k] = (*field_avg)[j][k] + grid_value(component,i,j,k); 
			}
		}

		if(t == end)
		{