input.c

program FDTD for the photonic crystal structure

}

float rand_distance_2D(float x_scan, float y_scan, float x_temp, float y_temp)
{
	return( sqrt((x_scan-x_temp)*(x_scan-x_temp) + (y_scan-y_temp)*(y_scan-y_temp)) );
}

void generate_random_3D(float x_min, float x_max, float y_min, float y_max, float z_min, float z_max, float radius, float *xr, float *yr, float *zr, int i)
{
	int scan; 
	float xr_temp, yr_temp, zr_temp; 	

	xr_temp = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(x_max-x_min) + x_min;	
	yr_temp = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(y_max-y_min) + y_min;
	zr_temp = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(z_max-z_min) + z_min;		

	if(i>0)
	{
		while(1)
		{
			for(scan=0; scan<i; scan++)
				if( rand_distance_3D(xr[scan], yr[scan], zr[scan], xr_temp, yr_temp, zr_temp) < 2*radius )
				{	
					scan = -10; 
					break;
				}

			if(scan == i)
			{
				xr[i] = xr_temp;
				yr[i] = yr_temp;
				zr[i] = zr_temp;
				break; //out of the while loop
			}
			else
			{
				xr_temp = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(x_max-x_min) + x_min;	
				yr_temp = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(y_max-y_min) + y_min;
				zr_temp = ( (float)rand() / ((float)(RAND_MAX)+(float)(1)) )*(z_max-z_min) + z_min;	
			}
		}
	}

	if(i==0)
	{
		xr[i] = xr_temp;
		yr[i] = yr_temp;
		zr[i] = zr_temp;
	}
}

float rand_distance_3D(float x_scan, float y_scan, float z_scan, float x_temp, float y_temp, float z_temp)
{
	return( sqrt((x_scan-x_temp)*(x_scan-x_temp) + (y_scan-y_temp)*(y_scan-y_temp) + (z_scan-z_temp)*(z_scan-z_temp)) );
}

void make_epsilon()
{
	int i,j,k;

	printf("-----------------------\n");	
	printf("making the epsilon...\n");

	for(i=0;i<misize;i++)
	{
		printf("%d% \n",100*i/misize); 

		for(j=0;j<mjsize;j++)
			for(k=0;k<mksize;k++)
			{
				epsilonx[i][j][k]=average_epsilon(i+0.5,j,k);
				epsilony[i][j][k]=average_epsilon(i,j+0.5,k);
				epsilonz[i][j][k]=average_epsilon(i,j,k+0.5);
			}
	}

	free(object);

	printf("make_epsilon...ok\n");
}

float average_epsilon(float i,float j,float k)
{
	int n,m,partial=0;
	float ii,jj,kk,epsilon,partial_epsilon[1000];

	epsilon=back_epsilon;
	for(m=0;m<1000;m++) partial_epsilon[m]=back_epsilon;

	for(n=0;n<objectn;n++)
	{
	
		if( in_object(n,i-0.5,j-0.5,k-0.5)==0 && in_object(n,i-0.5,j-0.5,k+0.5)==0 && 
			in_object(n,i-0.5,j+0.5,k-0.5)==0 && in_object(n,i-0.5,j+0.5,k+0.5)==0 &&
			in_object(n,i+0.5,j-0.5,k-0.5)==0 && in_object(n,i+0.5,j-0.5,k+0.5)==0 &&	
			in_object(n,i+0.5,j+0.5,k-0.5)==0 && in_object(n,i+0.5,j+0.5,k+0.5)==0);

		else if(in_object(n,i-0.5,j-0.5,k-0.5)==1 && in_object(n,i-0.5,j-0.5,k+0.5)==1 && 
				in_object(n,i-0.5,j+0.5,k-0.5)==1 && in_object(n,i-0.5,j+0.5,k+0.5)==1 &&
				in_object(n,i+0.5,j-0.5,k-0.5)==1 && in_object(n,i+0.5,j-0.5,k+0.5)==1 &&
				in_object(n,i+0.5,j+0.5,k-0.5)==1 && in_object(n,i+0.5,j+0.5,k+0.5)==1)
		{
			epsilon=object[n].epsilon;
			for(m=0;m<1000;m++) partial_epsilon[m]=object[n].epsilon;
			partial=0;
		}

		else
		{
			for(m=0,ii=i-0.45;ii<i+0.5;ii=ii+0.1)
			for(jj=j-0.45;jj<j+0.5;jj=jj+0.1)
			for(kk=k-0.45;kk<k+0.5;kk=kk+0.1,m++)
				if(in_object(n,ii,jj,kk)==1) partial_epsilon[m]=object[n].epsilon;
			partial=1;
		}
	}

	if(partial==1)
		for(epsilon=0,m=0;m<1000;m++)
			epsilon=epsilon+partial_epsilon[m]/1000;
	
	return epsilon;
}

int in_object(int n,float i,float j,float k)
{
	int I, J; 
	float X, Y, Z;

	if(strcmp(object[n].shape,"Rrod")!=0 && strcmp(object[n].shape,"Rellipse")!=0 && strcmp(object[n].shape,"Rellipsoidal")!=0 && strcmp(object[n].shape,"Rblock")!=0)
	{
		if(strcmp(object[n].shape,"rod")==0)
		{
			if( (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)<=object[n].size1*object[n].size1 && non_uniform_z_to_i(object[n].centerk-0.5*object[n].size2)<=k && k<=non_uniform_z_to_i(object[n].centerk+0.5*object[n].size2) ) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"rodX")==0)
		{
			if( (j-object[n].centerj)*(j-object[n].centerj)+(k-object[n].centerk)*(k-object[n].centerk)<=object[n].size1*object[n].size1 && object[n].centeri-object[n].size2/2<=i && i<=object[n].centeri+object[n].size2/2) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"rodY")==0)
		{
			if( (k-object[n].centerk)*(k-object[n].centerk)+(i-object[n].centeri)*(i-object[n].centeri)<=object[n].size1*object[n].size1 && object[n].centerj-object[n].size2/2<=j && j<=object[n].centerj+object[n].size2/2) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"donut")==0)
		{
			if( (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)>=object[n].size1*object[n].size1 && (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)<=object[n].size2*object[n].size2 && object[n].centerk-object[n].size3/2<=k && k<=object[n].centerk+object[n].size3/2) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"sphere")==0)
		{
			if( (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)+(k-object[n].centerk)*(k-object[n].centerk)<=object[n].size1*object[n].size1 ) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"shell")==0)
		{
			if( (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)+(k-object[n].centerk)*(k-object[n].centerk)>=object[n].size1*object[n].size1 && (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)+(k-object[n].centerk)*(k-object[n].centerk)<=object[n].size2*object[n].size2 ) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"ellipse")==0)
		{
			if( (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)/(object[n].size3*object[n].size3)<=object[n].size1*object[n].size1 && object[n].centerk-object[n].size2/2<=k && k<=object[n].centerk+object[n].size2/2) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"ellipsoidal")==0)
		{
			if( (i-object[n].centeri)*(i-object[n].centeri)/((object[n].size1)*(object[n].size1))+(j-object[n].centerj)*(j-object[n].centerj)/((object[n].size2)*(object[n].size2))+(k-object[n].centerk)*(k-object[n].centerk)/((object[n].size3)*(object[n].size3))<=1 ) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"cone")==0)
		{
			if( (i-object[n].centeri)*(i-object[n].centeri)+(j-object[n].centerj)*(j-object[n].centerj)<= ((object[n].size1-object[n].size2)*(k-object[n].centerk)/object[n].size3+0.5*(object[n].size1+object[n].size2))*((object[n].size1-object[n].size2)*(k-object[n].centerk)/object[n].size3+0.5*(object[n].size1+object[n].size2)) && (k-object[n].centerk)<= 0.5*object[n].size3 && (k-object[n].centerk)>= -0.5*object[n].size3) return 1;
			else return 0;
		}

		else if(strcmp(object[n].shape,"block")==0)
		{
			if( object[n].centeri-object[n].size1/2<=i && i<=object[n].centeri+object[n].size1/2 &&
				object[n].centerj-object[n].size2/2<=j && j<=object[n].centerj+object[n].size2/2 &&
				non_uniform_z_to_i(object[n].centerk-0.5*object[n].size3)<=k && k<=non_uniform_z_to_i(object[n].centerk+0.5*object[n].size3) ) return 1;
			else return 0;
		}

		else if(strcmp(object[n].shape,"contour")==0)
		{
			if( (-object[n].centeri*CUTFACTOR <= (i-(xcenter*lattice_x))/object[n].size3) && ((i-(xcenter*lattice_x))/object[n].size3 <= (object[n].col-object[n].centeri)*CUTFACTOR) && (-object[n].centerj*CUTFACTOR <= (j-(ycenter*lattice_y))/object[n].size3) && ((j-(ycenter*lattice_y))/object[n].size3 <= (object[n].row-object[n].centerj)*CUTFACTOR) )
			{
				I = floor( object[n].centeri + (i - (xcenter*lattice_x))/object[n].size3 ); 
				J = floor( object[n].centerj + (j - (ycenter*lattice_y))/object[n].size3 );
		       		if( (object[n].matrix[I][J] <= object[n].size1) && (object[n].centerk-object[n].size2/2 <= k) && (k <= object[n].centerk+object[n].size2/2) ) 
				return 1;
			}
			else return 0;
		}
	}
	else 	
////if(strcmp(object[n].shape,"Rrod")==0 || strcmp(object[n].shape,"Rellipse")==0 || 
////strcmp(object[n].shape,"Rellipsoidal")==0 || strcmp(object[n].shape,"Rblock")==0)
	{
		//// Transform (i,j,k) --> (X,Y,Z), the coordinate with respect to the object
		X = (cos_a[n]*cos_c[n]-sin_a[n]*cos_b[n]*sin_c[n])*(i-object[n].centeri) 
			+ (sin_a[n]*cos_c[n]+cos_a[n]*cos_b[n]*sin_c[n])*(j-object[n].centerj)
			+ (sin_b[n]*sin_c[n])*(k-object[n].centerk) + object[n].centeri;
		Y =  -(cos_a[n]*sin_c[n]+sin_a[n]*cos_b[n]*cos_c[n])*(i-object[n].centeri)
			-(sin_a[n]*sin_c[n]-cos_a[n]*cos_b[n]*cos_c[n])*(j-object[n].centerj)
			+(sin_b[n]*cos_c[n])*(k-object[n].centerk) + object[n].centerj;
		Z = (sin_a[n]*sin_b[n])*(i-object[n].centeri) - (cos_a[n]*sin_b[n])*(j-object[n].centerj)
		 	+ cos_b[n]*(k-object[n].centerk) + object[n].centerk;

		if(strcmp(object[n].shape,"Rrod")==0)
		{
			if( (X-object[n].centeri)*(X-object[n].centeri)+(Y-object[n].centerj)*(Y-object[n].centerj)<=object[n].size1*object[n].size1 && object[n].centerk-object[n].size2/2<=Z && Z<=object[n].centerk+object[n].size2/2) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"Rellipse")==0)
		{
			if( (X-object[n].centeri)*(X-object[n].centeri)+(Y-object[n].centerj)*(Y-object[n].centerj)/(object[n].size3*object[n].size3)<=object[n].size1*object[n].size1 && object[n].centerk-object[n].size2/2<=Z && Z<=object[n].centerk+object[n].size2/2) return 1;
			else return 0;
		}

		if(strcmp(object[n].shape,"Rellipsoidal")==0)
		{
			if( (X-object[n].centeri)*(X-object[n].centeri)/((object[n].size1)*(object[n].size1))+(Y-object[n].centerj)*(Y-object[n].centerj)/((object[n].size2)*(object[n].size2))+(Z-object[n].centerk)*(Z-object[n].centerk)/((object[n].size3)*(object[n].size3))<=1 ) return 1;
			else return 0;
		}

		else if(strcmp(object[n].shape,"Rblock")==0)
		{
			if( object[n].centeri-object[n].size1/2<=X && X<=object[n].centeri+object[n].size1/2 &&
				object[n].centerj-object[n].size2/2<=Y && Y<=object[n].centerj+object[n].size2/2 &&
				object[n].centerk-object[n].size3/2<=Z && Z<=object[n].centerk+object[n].size3/2) return 1;
			else return 0;
		}
	}
}		   

void reading_matrix(char *matrix_file, int n)
{
	FILE *stream;
	int i,j;
	int col, row;
	char ch;
	char string[20];
	
	stream = fopen(matrix_file,"rt");

	col=0; row=0;	
	/////// Counting row and col ///////
	ch = getc(stream);
	while( ch != EOF )
	{
		if( ch == '\t' ) 
			col++;
		if( ch == '\n' )
		{	
			col++;    // Origin standard (with TAB separation)
			row++;
		}
		ch = getc(stream);
	}
	col = (int)(col / row); // in each row
	printf("%s matrix file\n",matrix_file);
	printf("  matrix col = %d\n", col);
	printf("  matrix row = %d\n", row);
	printf("\n");
	fclose(stream);

	object[n].col=col; object[n].row=row;

	//////// matrix memory allocation /////
	object[n].matrix = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		object[n].matrix[i] = (float *)malloc(sizeof(float)*row);
	
	stream = fopen(matrix_file,"rt");
	//////// Reading matrix file data //////
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream,"%s",string);
			object[n].matrix[i][j] = atof(string);
		}
	}

	fclose(stream);	

}