output.c

program FDTD for the photonic crystal structure

{
	int deltar; // (NROW-row)/2
	int deltac;
	char ch;
	char string[80];  

	double **source; // source matrix

	int row, col; // row num & column num
	int i,j;

	FILE *stream;
	char name_freq[10];

	char name_tEx_real[20], name_tEx_imag[20], name_tEy_real[20], name_tEy_imag[20];
	char name_tHx_real[20], name_tHx_imag[20], name_tHy_real[20], name_tHy_imag[20];

	char name_Ex_real[20], name_Ex_imag[20], name_Ey_real[20], name_Ey_imag[20];
	char name_Hx_real[20], name_Hx_imag[20], name_Hy_real[20], name_Hy_imag[20];

	/////////////////////////////////////////////////
	//// Making file names   ////////////////////////
	/////////////////////////////////////////////////
	sprintf(name_freq,".ri%02d",mm);

	//// Source file name ////
	////----------------------
	sprintf(name_tEx_real,"tEx_real");
	sprintf(name_tEx_imag,"tEx_imag");
	strcat(name_tEx_real,name_freq);
	strcat(name_tEx_imag,name_freq);

	sprintf(name_tEy_real,"tEy_real");
	sprintf(name_tEy_imag,"tEy_imag");
	strcat(name_tEy_real,name_freq);
	strcat(name_tEy_imag,name_freq);

	sprintf(name_tHx_real,"tHx_real");
	sprintf(name_tHx_imag,"tHx_imag");
	strcat(name_tHx_real,name_freq);
	strcat(name_tHx_imag,name_freq);

	sprintf(name_tHy_real,"tHy_real");
	sprintf(name_tHy_imag,"tHy_imag");
	strcat(name_tHy_real,name_freq);
	strcat(name_tHy_imag,name_freq);

	//// Output file name ////
	////----------------------
	sprintf(name_Ex_real,"Ex_real");
	sprintf(name_Ex_imag,"Ex_imag");
	strcat(name_Ex_real,name_freq);
	strcat(name_Ex_imag,name_freq);

	sprintf(name_Ey_real,"Ey_real");
	sprintf(name_Ey_imag,"Ey_imag");
	strcat(name_Ey_real,name_freq);
	strcat(name_Ey_imag,name_freq);

	sprintf(name_Hx_real,"Hx_real");
	sprintf(name_Hx_imag,"Hx_imag");
	strcat(name_Hx_real,name_freq);
	strcat(name_Hx_imag,name_freq);

	sprintf(name_Hy_real,"Hy_real");
	sprintf(name_Hy_imag,"Hy_imag");
	strcat(name_Hy_real,name_freq);
	strcat(name_Hy_imag,name_freq);

	/////////////////////////////////////////////////
	// Count row and col ////////////////////////////
	/////////////////////////////////////////////////

	col = row = 0; 
	stream = fopen(name_tEx_real,"r");
	ch = getc(stream);
	while( ch != EOF )
	{
		if( ch == '\t' ) 
			col++;
		if( ch == '\n' )
		{
			row++;
		}
		ch = getc(stream);
	}
	fclose(stream);
	col = (int)(col / row); // in each row
	printf("t_data col = %d\n", col);
	printf("t_data row = %d\n", row);

	/////////////////////////////////////////////////
	// reading source ///////////////////////////////
	/////////////////////////////////////////////////

	deltar = (NROW-row)/2;
	deltac = (NROW-col)/2;

	source = (double **)malloc(sizeof(double *)*col);
	for(i=0; i<col; i++)
		source[i] = (double *)malloc(sizeof(double)*row);

	stream = fopen(name_tEx_real,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Ex_real,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	stream = fopen(name_tEx_imag,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Ex_imag,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	stream = fopen(name_tEy_real,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Ey_real,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	stream = fopen(name_tEy_imag,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Ey_imag,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	///////////////////////////////////////////

	stream = fopen(name_tHx_real,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Hx_real,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	stream = fopen(name_tHx_imag,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Hx_imag,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	stream = fopen(name_tHy_real,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Hy_real,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	stream = fopen(name_tHy_imag,"r");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", &string);
			source[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Hy_imag,"wt");
	for(j=NROW-1; j>=0; j--)
	{
		for(i=0; i<NROW; i++)
		{
			if( (i>deltac && i<(NROW-deltac)) && (j>deltar && j<(NROW-deltar)) )
				fprintf(stream, "%lf\t", source[i-deltac][j-deltar]);
			else if( i == (NROW-1) )
				fprintf(stream, "%lf\t\n", 0.0);
			else
				fprintf(stream, "%lf\t", 0.0);				
		}
	}
	fclose(stream);

	free(source);
}

void print_real_and_imag(int mm)
{
	int i, j;
	FILE *real, *imag;
	char name_freq[10];

	char name_Ex_real[20], name_Ex_imag[20], name_Ey_real[20], name_Ey_imag[20];
	char name_Hx_real[20], name_Hx_imag[20], name_Hy_real[20], name_Hy_imag[20];

	//////// making file names /////////
	sprintf(name_freq,".ri%02d",mm);

	sprintf(name_Ex_real,"tEx_real");
	sprintf(name_Ex_imag,"tEx_imag");
	strcat(name_Ex_real,name_freq);
	strcat(name_Ex_imag,name_freq);

	sprintf(name_Ey_real,"tEy_real");
	sprintf(name_Ey_imag,"tEy_imag");
	strcat(name_Ey_real,name_freq);
	strcat(name_Ey_imag,name_freq);

	sprintf(name_Hx_real,"tHx_real");
	sprintf(name_Hx_imag,"tHx_imag");
	strcat(name_Hx_real,name_freq);
	strcat(name_Hx_imag,name_freq);

	sprintf(name_Hy_real,"tHy_real");
	sprintf(name_Hy_imag,"tHy_imag");
	strcat(name_Hy_real,name_freq);
	strcat(name_Hy_imag,name_freq);

	/////////////////////////////////////
	/////////// Print Ex dat ////////////
	/////////////////////////////////////
	real = fopen(name_Ex_real,"wt");
	imag = fopen(name_Ex_imag,"wt");
	for(j=jsize-2; j>=1; j--)
	{
		for(i=1; i<=isize-2; i++)
		{
			if( (i >= pmlil && i <=isize-pmlir) && (j >= pmljl && j <= jsize-pmljr))
			{
				fprintf(real,"%f\t", Ex_cos[i][j][mm]);
				fprintf(imag,"%f\t", Ex_sin[i][j][mm]);
			}
			else
			{
				fprintf(real,"%f\t", 0.0);
				fprintf(imag,"%f\t", 0.0);
			}
		}
		fprintf(real, "\n");
		fprintf(imag, "\n");
	}
	fclose(real); fclose(imag);

	/////////////////////////////////////
	/////////// Print Ey dat ////////////
	/////////////////////////////////////
	real = fopen(name_Ey_real,"wt");
	imag = fopen(name_Ey_imag,"wt");
	for(j=jsize-2; j>=1; j--)
	{
		for(i=1; i<=isize-2; i++)
		{
			if( (i >= pmlil && i <=isize-pmlir) && (j >= pmljl && j <= jsize-pmljr))
			{
				fprintf(real,"%f\t", Ey_cos[i][j][mm]);
				fprintf(imag,"%f\t", Ey_sin[i][j][mm]);
			}
			else
			{
				fprintf(real,"%f\t", 0.0);
				fprintf(imag,"%f\t", 0.0);
			}
		}
		fprintf(real, "\n");
		fprintf(imag, "\n");
	}
	fclose(real); fclose(imag);

	/////////////////////////////////////
	/////////// Print Hx dat ////////////
	/////////////////////////////////////
	real = fopen(name_Hx_real,"wt");
	imag = fopen(name_Hx_imag,"wt");
	for(j=jsize-2; j>=1; j--)
	{
		for(i=1; i<=isize-2; i++)
		{
			if( (i >= pmlil && i <=isize-pmlir) && (j >= pmljl && j <= jsize-pmljr))
			{
				fprintf(real,"%f\t", Hx_cos[i][j][mm]);
				fprintf(imag,"%f\t", Hx_sin[i][j][mm]);
			}
			else
			{
				fprintf(real,"%f\t", 0.0);
				fprintf(imag,"%f\t", 0.0);
			}
		}
		fprintf(real, "\n");
		fprintf(imag, "\n");
	}
	fclose(real); fclose(imag);

	/////////////////////////////////////
	/////////// Print Hy dat ////////////
	/////////////////////////////////////
	real = fopen(name_Hy_real,"wt");
	imag = fopen(name_Hy_imag,"wt");
	for(j=jsize-2; j>=1; j--)
	{
		for(i=1; i<=isize-2; i++)
		{
			if( (i >= pmlil && i <=isize-pmlir) && (j >= pmljl && j <= jsize-pmljr))
			{
				fprintf(real,"%f\t", Hy_cos[i][j][mm]);
				fprintf(imag,"%f\t", Hy_sin[i][j][mm]);
			}
			else
			{
				fprintf(real,"%f\t", 0.0);
				fprintf(imag,"%f\t", 0.0);
			}
		}
		fprintf(real, "\n");
		fprintf(imag, "\n");
	}
	fclose(real); fclose(imag);

	printf("print real and imag [%d] ok....!\n",mm);

}

int find_k_projection(int i, int j, int k_range)
{
	float initial;
	int m;

	initial = eps(i,j,k_range);

	for(m=k_range-1; m>0; m--)
	{
		if( meps(i,j,m)==0.0 && fabs(eps(i,j,m)-initial)>0.1 )
			return(m-1);
		if( meps(i,j,m)!=0.0 )
			return(m-1);
	}
	return(0);
}