output.c

program FDTD for the photonic crystal structure

	return p*(Jy[i][j][k]+Jy[i][j-1][k])/2;
}

float jfieldz(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(Jz[i][j][k]+Jz[i][j][k-1])/2;
}

float hfieldx(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(Hx[i][j][k]+Hx[i][j-1][k]+Hx[i][j][k-1]+Hx[i][j-1][k-1])/4;
}

float hfieldy(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(Hy[i][j][k]+Hy[i][j][k-1]+Hy[i-1][j][k]+Hy[i-1][j][k-1])/4;
}

float hfieldz(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(Hz[i][j][k]+Hz[i-1][j][k]+Hz[i][j-1][k]+Hz[i-1][j-1][k])/4;
}

float iefieldx(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(iEx[i][j][k]+iEx[i-1][j][k])/2;
}

float iefieldy(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(iEy[i][j][k]+iEy[i][j-1][k])/2;
}

float iefieldz(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(iEz[i][j][k]+iEz[i][j][k-1])/2;
}

float ijfieldx(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(iJx[i][j][k]+iJx[i-1][j][k])/2;
}

float ijfieldy(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(iJy[i][j][k]+iJy[i][j-1][k])/2;
}

float ijfieldz(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(iJz[i][j][k]+iJz[i][j][k-1])/2;
}

float ihfieldx(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(iHx[i][j][k]+iHx[i][j-1][k]+iHx[i][j][k-1]+iHx[i][j-1][k-1])/4;
}

float ihfieldy(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(iHy[i][j][k]+iHy[i][j][k-1]+iHy[i-1][j][k]+iHy[i-1][j][k-1])/4;
}

float ihfieldz(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(iHz[i][j][k]+iHz[i-1][j][k]+iHz[i][j-1][k]+iHz[i-1][j-1][k])/4;
}

float eps(int i,int j,int k)
{
	if(pisize<=i && (xparity==1 || xparity==-1)) i=isize-i;  
	if(pjsize<=j && (yparity==1 || yparity==-1)) j=jsize-j; 
	if(pksize<=k && (zparity==1 || zparity==-1)) k=ksize-k;

	return (epsilonx[i][j][k]+epsilonx[i-1][j][k]+epsilony[i][j][k]+epsilony[i][j-1][k]+epsilonz[i][j][k]+epsilonz[i][j][k-1])/6;
}

float meps(int i,int j,int k)
{
	if(pisize<=i && (xparity==1 || xparity==-1)) i=isize-i;  
	if(pjsize<=j && (yparity==1 || yparity==-1)) j=jsize-j; 
	if(pksize<=k && (zparity==1 || zparity==-1)) k=ksize-k;

	return mepsilon[i][j][k];
}

void print_amp_and_phase(int mm)
{
	int i, j;
	FILE *amp, *phase;
	char name_freq[10];

	char name_Ex_amp[20], name_Ex_phase[20], name_Ey_amp[20], name_Ey_phase[20];
	char name_Hx_amp[20], name_Hx_phase[20], name_Hy_amp[20], name_Hy_phase[20];

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

	sprintf(name_Ex_amp,"jEx_amp");
	sprintf(name_Ex_phase,"jEx_phase");
	strcat(name_Ex_amp,name_freq);
	strcat(name_Ex_phase,name_freq);

	sprintf(name_Ey_amp,"jEy_amp");
	sprintf(name_Ey_phase,"jEy_phase");
	strcat(name_Ey_amp,name_freq);
	strcat(name_Ey_phase,name_freq);

	sprintf(name_Hx_amp,"jHx_amp");
	sprintf(name_Hx_phase,"jHx_phase");
	strcat(name_Hx_amp,name_freq);
	strcat(name_Hx_phase,name_freq);

	sprintf(name_Hy_amp,"jHy_amp");
	sprintf(name_Hy_phase,"jHy_phase");
	strcat(name_Hy_amp,name_freq);
	strcat(name_Hy_phase,name_freq);

	/////////////////////////////////////
	/////////// Print Ex dat ////////////
	/////////////////////////////////////
	amp = fopen(name_Ex_amp,"wt");
	phase = fopen(name_Ex_phase,"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(amp,"%f\t", sqrt( pow(Ex_cos[i][j][mm],2) + pow(Ex_sin[i][j][mm],2) ));
				fprintf(phase,"%f\t", -atan2( Ex_sin[i][j][mm], Ex_cos[i][j][mm] ));
			}
			else
			{
				fprintf(amp,"%f\t", 0.0);
				fprintf(phase,"%f\t", 0.0);
			}
		}
		fprintf(amp, "\n");
		fprintf(phase, "\n");
	}
	fclose(amp); fclose(phase);

	/////////////////////////////////////
	/////////// Print Ey dat ////////////
	/////////////////////////////////////
	amp = fopen(name_Ey_amp,"wt");
	phase = fopen(name_Ey_phase,"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(amp,"%f\t", sqrt( pow(Ey_cos[i][j][mm],2) + pow(Ey_sin[i][j][mm],2) ));
				fprintf(phase,"%f\t", -atan2( Ey_sin[i][j][mm], Ey_cos[i][j][mm] ));
			}
			else
			{
				fprintf(amp, "%f\t", 0.0);
				fprintf(phase, "%f\t", 0.0);
			}
		}
		fprintf(amp, "\n");
		fprintf(phase, "\n");
	}
	fclose(amp); fclose(phase);

	/////////////////////////////////////
	/////////// Print Hx dat ////////////
	/////////////////////////////////////
	amp = fopen(name_Hx_amp,"wt");
	phase = fopen(name_Hx_phase,"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(amp,"%f\t", sqrt( pow(Hx_cos[i][j][mm],2) + pow(Hx_sin[i][j][mm],2) ));
				fprintf(phase,"%f\t", -atan2( Hx_sin[i][j][mm], Hx_cos[i][j][mm] ));
			}
			else
			{
				fprintf(amp, "%f\t", 0.0);
				fprintf(phase, "%f\t", 0.0);
			}
		}
		fprintf(amp, "\n");
		fprintf(phase, "\n");
	}
	fclose(amp); fclose(phase);

	/////////////////////////////////////
	/////////// Print Hy dat ////////////
	/////////////////////////////////////
	amp = fopen(name_Hy_amp,"wt");
	phase = fopen(name_Hy_phase,"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(amp,"%f\t", sqrt( pow(Hy_cos[i][j][mm],2) + pow(Hy_sin[i][j][mm],2) ));
				fprintf(phase,"%f\t", -atan2( Hy_sin[i][j][mm], Hy_cos[i][j][mm] ));
			}
			else
			{
				fprintf(amp, "%f\t", 0.0);
				fprintf(phase, "%f\t", 0.0);
			}
		}
		fprintf(amp, "\n");
		fprintf(phase, "\n");
	}
	fclose(amp); fclose(phase);

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

}

void print_real_and_imag_2n_size(int NROW, int mm)
{
	char string[80];  

	int row, col; // row num & column num
	int deltar; // (nrow-row)/2
	int deltac;

	int i,j;
	FILE *stream;

	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,"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);

	/////////////////////////////////////////////////
	// counting row num, col num ////////////////////
	/////////////////////////////////////////////////
	row = jsize-2; col = isize-2; 
	deltar = (NROW-row)/2 + (pmljl+10+pmljr+10)/2;  // triming out pml region that may cause numerical errors
	deltac = (NROW-col)/2 + (pmlil+10+pmlir+10)/2;

	/////////////////////////////////////////////////
	// File out enlarged filed data /////////////////
	/////////////////////////////////////////////////

	///////// Ex component ////////////	
	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, "%f\t", Ex_cos[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	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, "%f\t", Ex_sin[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	fclose(stream);

	///////// Ey component ////////////	
	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, "%f\t", Ey_cos[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	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, "%f\t", Ey_sin[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	fclose(stream);

	///////// Hx component ////////////	
	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, "%f\t", Hx_cos[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	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, "%f\t", Hx_sin[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	fclose(stream);

	///////// Hy component ////////////	
	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, "%f\t", Hy_cos[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	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, "%f\t", Hy_sin[i-deltac+(pmlil+10)][j-deltar+(pmlir+10)][mm] );
			else
				fprintf(stream, "%f\t", 0.0);				
		}
		fprintf(stream, "\n");
	}
	fclose(stream);

	printf("make 2n_size data [%d] ok....!\n",mm);
}

void make_2n_size(int NROW, int mm)