input.c

program FDTD for the photonic crystal structure

#include "./pFDTD.h"

int in_object(int n,float i,float j,float k);
void reading_matrix(char *matrix_file, int n);
float average_epsilon(float i,float j,float k);
void generate_random_2D(float x_min, float x_max, float y_min, float y_max, float radius, float *xr, float *yr, int i);
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);
float rand_distance_2D(float x_scan, float y_scan, float x_temp, float y_temp);
float rand_distance_3D(float x_scan, float y_scan, float z_scan, float x_temp, float y_temp, float z_temp);

float back_epsilon;
int objectn=0;
struct obj *object;

///// For Euler rotation ///////
///// [n] : object index ///////
/// a = alpha, b=beta, c=gamma ///
float *cos_a, *cos_b, *cos_c;
float *sin_a, *sin_b, *sin_c;
//////////////////////////////////////

void background(float epsilon)
{
	int i,j,k;

	back_epsilon=epsilon;
	
	for(i=0;i<misize;i++)
	for(j=0;j<mjsize;j++)
	for(k=0;k<mksize;k++)
		epsilonx[i][j][k]=epsilony[i][j][k]=epsilonz[i][j][k]=back_epsilon;

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

void input_object(char *shape, char *matrix_file, float centerx, float centery, float centerz, float size1, float size2, float size3, float epsilon)
{	
	objectn++;

	object=(struct obj *)realloc(object,sizeof(struct obj)*objectn);
	
	cos_a=(float *)realloc(cos_a,sizeof(float)*objectn);
	cos_b=(float *)realloc(cos_b,sizeof(float)*objectn);
	cos_c=(float *)realloc(cos_c,sizeof(float)*objectn);
	sin_a=(float *)realloc(sin_a,sizeof(float)*objectn);
	sin_b=(float *)realloc(sin_b,sizeof(float)*objectn);
	sin_c=(float *)realloc(sin_c,sizeof(float)*objectn);
	
	strcpy(object[objectn-1].shape,shape);

	//////// common parameters ///////
	object[objectn-1].epsilon=epsilon;  
	//////////////////////////////////
	cos_a[objectn-1] = 1.0;
	cos_b[objectn-1] = 1.0;
	cos_c[objectn-1] = 1.0;
	sin_a[objectn-1] = 0.0;
	sin_b[objectn-1] = 0.0;
	sin_c[objectn-1] = 0.0;

	// in the below, we define various input parameters case by case. 

	if(strcmp(shape,"contour")==0)
	{
		object[objectn-1].centeri=centerx;           // matrix_center 
		object[objectn-1].centerj=centery;           // matrix_center 
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1;     // discrimination level
		object[objectn-1].size2=size2*lattice_z;   // height
		object[objectn-1].size3=size3;     // compression factor
		reading_matrix(matrix_file, objectn-1);	
	}
	else if(strcmp(shape,"rod")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;  
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz); // for non_uniform_grid(), the meaning changed
		object[objectn-1].size1=size1*lattice_x;  // radius
		object[objectn-1].size2=(size2);  // height  // for non_uniform_grid() ....
	}
	else if(strcmp(shape,"rodX")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;  
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_y;  // radius
		object[objectn-1].size2=size2*lattice_x;  // x length
	}
	else if(strcmp(shape,"rodY")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;  
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // radius
		object[objectn-1].size2=size2*lattice_y;  // y length
	}
	else if(strcmp(shape,"donut")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;  
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // inner radius
		object[objectn-1].size2=size2*lattice_x;  // outer radius
		object[objectn-1].size3=size3*lattice_z;  // height
	}
	else if(strcmp(shape,"sphere")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // radius
	}
	else if(strcmp(shape,"shell")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // inner radius
		object[objectn-1].size2=size2*lattice_x;  // outer radius
	}
	else if(strcmp(shape,"ellipse")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // radius
		object[objectn-1].size2=size2*lattice_z;  // height
		object[objectn-1].size3=size3;  // aspect ratio=ry/rx
	}
	else if(strcmp(shape,"ellipsoidal")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // rx
		object[objectn-1].size2=size2*lattice_y;  // ry
		object[objectn-1].size3=size3*lattice_z;  // rz
	}
	else if(strcmp(shape,"cone")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // r1
		object[objectn-1].size2=size2*lattice_x;  // r2
		object[objectn-1].size3=size3*lattice_z;  // height
	}
	else // "block"
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz); // for non_uniform_grid(), the meaning changed
		object[objectn-1].size1=size1*lattice_x;
		object[objectn-1].size2=size2*lattice_y;
		object[objectn-1].size3=(size3);  // for non_uniform_grid(), the meaning changed
	}
}

void input_object_Euler_rotation(char *shape, char *matrix_file, float centerx, float centery, float centerz, float size1, float size2, float size3, float alpha, float beta, float gamma, float epsilon)
// alpha, beta, gamma : 3- Euler angles 
{	
	objectn++;

	object=(struct obj *)realloc(object,sizeof(struct obj)*objectn);

	cos_a=(float *)realloc(cos_a,sizeof(float)*objectn);
	cos_b=(float *)realloc(cos_b,sizeof(float)*objectn);
	cos_c=(float *)realloc(cos_c,sizeof(float)*objectn);
	sin_a=(float *)realloc(sin_a,sizeof(float)*objectn);
	sin_b=(float *)realloc(sin_b,sizeof(float)*objectn);
	sin_c=(float *)realloc(sin_c,sizeof(float)*objectn);
	
	strcpy(object[objectn-1].shape,shape);

	//////// common parameters ///////
	object[objectn-1].epsilon=epsilon;  
	//////////////////////////////////
	cos_a[objectn-1] = cos(alpha*pi/180);
	cos_b[objectn-1] = cos(beta*pi/180);
	cos_c[objectn-1] = cos(gamma*pi/180);
	sin_a[objectn-1] = sin(alpha*pi/180);
	sin_b[objectn-1] = sin(beta*pi/180);
	sin_c[objectn-1] = sin(gamma*pi/180);

	//// in the below, we define various input parameters case by case. ////
	//// Contrary to "input_object()", 'Euler rotaion' cannot be applied to "sphere" and "contour"////

	if(strcmp(shape,"rod")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;  
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // radius
		object[objectn-1].size2=size2*lattice_z;  // height
		strcpy(object[objectn-1].shape,"Rrod"); // Euler Rotation indicator 
	}
	else if(strcmp(shape,"ellipse")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // radius
		object[objectn-1].size2=size2*lattice_z;  // height
		object[objectn-1].size3=size3;  // aspect ratio=ry/rx
		strcpy(object[objectn-1].shape,"Rellipse"); // Euler Rotation indicator 
	}
	else if(strcmp(shape,"ellipsoidal")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;  // rx
		object[objectn-1].size2=size2*lattice_y;  // ry
		object[objectn-1].size3=size3*lattice_z;  // rz
		strcpy(object[objectn-1].shape,"Rellipsoidal"); // Euler Rotation indicator 
	}
	else if(strcmp(shape,"block")==0)
	{
		object[objectn-1].centeri=(centerx+xcenter)*lattice_x;
		object[objectn-1].centerj=(centery+ycenter)*lattice_y;
		object[objectn-1].centerk=(centerz+zcenter)*lattice_z;
		object[objectn-1].size1=size1*lattice_x;
		object[objectn-1].size2=size2*lattice_y;
		object[objectn-1].size3=size3*lattice_z;
		strcpy(object[objectn-1].shape,"Rblock"); // Euler Rotation indicator 
	}
	else
	{}
}

void random_object(char *shape, float radius, float height, float epsilon, float x_min, float x_max, float y_min, float y_max, float z_min, float z_max, int gen_number, int seed)
{
	int i;	
	float *xr, *yr, *zr;  //random number coordinates

	srand(seed); //initialization of random number generator

	if(strcmp(shape,"rod")==0)
	{
		xr = (float *)malloc(sizeof(float)*gen_number);
		yr = (float *)malloc(sizeof(float)*gen_number);

		for(i=0; i<gen_number; i++)
		{
			generate_random_2D(x_min, x_max, y_min, y_max, radius, xr, yr, i);
			input_object("rod", EMP, xr[i], yr[i], 0+shift, radius, height, 0, epsilon);
		}
	}

	if(strcmp(shape,"sphere")==0)
	{
		xr = (float *)malloc(sizeof(float)*gen_number);
		yr = (float *)malloc(sizeof(float)*gen_number);
		zr = (float *)malloc(sizeof(float)*gen_number);

		for(i=0; i<gen_number; i++)
		{
			generate_random_3D(x_min, x_max, y_min, y_max, z_min, z_max, radius, xr, yr, zr, i);
			input_object("sphere", EMP, xr[i], yr[i], zr[i], radius, 0, 0, epsilon);
		}
	}	
}

void generate_random_2D(float x_min, float x_max, float y_min, float y_max, float radius, float *xr, float *yr, int i)
{
	int scan; 
	float xr_temp, yr_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;	

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

			if(scan == i)
			{
				xr[i] = xr_temp;
				yr[i] = yr_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;	
			}
		}
	}

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