parameter.c

program FDTD for the photonic crystal structure

#include "./pFDTD.h"

int isize, jsize, ksize;
int pmlil, pmlir, pmljl, pmljr, pmlkl, pmlkr;
int lattice_x, lattice_y, lattice_z;
int lattice_nz;
float xsize, ysize, zsize;
float nz_start, nz_end;
float xcenter, ycenter, zcenter;
float kx, ky, kz;
float orderxl, orderyl, orderzl;
float orderxr, orderyr, orderzr;
float sig_axl, sig_ayl, sig_azl;
float sig_axr, sig_ayr, sig_azr;
float ds_x, ds_y, ds_z, dt;  
float *ds_nz;
float S_factor; 
float pi, eo, uo, ups, light_speed;
int misize, mjsize, mksize;
int pisize, pjsize, pksize;
int cisize, cjsize, cksize;
int xparity=0, yparity=0, zparity=0;  //default 

float wave_vector_x, wave_vector_y;
int use_periodic_x=0, use_periodic_y=0; //default  

int lcm3(int a, int b, int c);
int lcm2(int m, int n);
int gcd2(int m, int n);
float *float_1d_memory_for_parameter(int imax);

void structure_size(float x,float y,float z)
{
	xsize=x;
	ysize=y;
	zsize=z;
}

void lattice_size(int lx, int ly, int lz)
{
	lattice_x=lx;
	lattice_y=ly;
	lattice_z=lz;

	///////////////////////////////////////////
	//// by default for non_uniform_grid() ////
	///////////////////////////////////////////
	lattice_nz=lattice_z;
	nz_start = -1234; // some large negative value outside the range
	nz_end = -1234;
}

void non_uniform_grid(char *component, float z_i, float z_f, int nlz)
{
	// only for z direction
	lattice_nz=nlz;
	nz_start = z_i;
	nz_end = z_f; 
}

int non_uniform_z_to_i(float z)
{
	if( z <= nz_start ) 
		return( floor(0.5+(z+0.5*zsize)*lattice_z) );
	else if( (z > nz_start) && (z <= nz_end) )
		return( floor(0.5+( (nz_start+0.5*zsize)*lattice_z + (z-nz_start)*lattice_nz )) );
	else // (z > nz_end)
		return( floor(0.5+( (nz_start+0.5*zsize)*lattice_z + (nz_end-nz_start)*lattice_nz + (z-nz_end)*lattice_z )) );
}

void pml_size(int il,int ir,int jl,int jr,int kl,int kr)
{
	pmlil=il;
	pmlir=ir;
	pmljl=jl;
	pmljr=jr;
	pmlkl=kl;
	pmlkr=kr;
	//See the range of position[i][j][k] in memory.c.  
	if(il==0)
		pmlil=0;
	if(ir==0)
		pmlir=0;
	if(jl==0)
		pmljl=0;
	if(jr==0)
		pmljr=0;
	if(kl==0)
		pmlkl=0;
	if(kr==0)
		pmlkr=0;
}

void set_default_parameter(float S)
{
	int lcm_temp;
	int i;

	pi=3.141592;
	eo=8.854e-12;
	uo=pi*4.0e-7;
	ups=1.0;
	light_speed=1.0/sqrt(eo*uo);

	lcm_temp = lcm3(lattice_x, lattice_y, lattice_nz);

	ds_x = lcm_temp/lattice_x;
	ds_y = lcm_temp/lattice_y;
	///////////////////////////////////////////////////////////////////////////
	//in case of using non_uniform_grid(), ds_z takes the smallest feature size
	/////////////////////////////////////////////////////////////////////////// 
	ds_z = lcm_temp/lattice_nz;
	
	S_factor = S;

	dt=1/light_speed/S_factor; // S=Couriant parameter for numerical stability

	orderxl = 3.5; orderyl = 3.5; orderzl = 3.5;
	orderxr = 3.5; orderyr = 3.5; orderzr = 3.5;

	sig_axl = 1.0; sig_ayl = 1.0; sig_azl = 1.0;
	sig_axr = 1.0; sig_ayr = 1.0; sig_azr = 1.0;

	kx = 1.0; ky = 1.0; kz =1.0;
	
	isize=floor(0.5+(xsize*lattice_x));
	jsize=floor(0.5+(ysize*lattice_y));
	////////////////////////////////////
	///// for non_uniform_grid() ///////
	////////////////////////////////////
	ksize = non_uniform_z_to_i(0.5*zsize);
	// above expression is equivalent as ...
	//ksize=floor(0.5+((nz_start+0.5*zsize)*lattice_z + (nz_end-nz_start)*lattice_nz + (0.5*zsize-nz_end)*lattice_z));

	xcenter=xsize/2;
	ycenter=ysize/2;
	zcenter=zsize/2;  // always true even with non_uniform_grid()

	misize=isize;
	pisize=isize-1;
	cisize=isize/2;
	mjsize=jsize;
	pjsize=jsize-1;
	cjsize=jsize/2;
	mksize=ksize;
	pksize=ksize-1;
	cksize=ksize/2;  // in case of using non_uniform_grid(), this will not be used. 

	///////////////////////////////////////////////////
	//// create non-uniform z-grid map (ds_nz[i]) /////
	///////////////////////////////////////////////////
	ds_nz = float_1d_memory_for_parameter(mksize);
	for(i=0; i<mksize; i++)
		ds_nz[i] = lcm_temp/lattice_z;
	for(i=non_uniform_z_to_i(nz_start); i<non_uniform_z_to_i(nz_end); i++)
		ds_nz[i] = lcm_temp/lattice_nz;
	///////////////////////////////////////////////////

	////// Random number generation //////
	gsl_rng_env_setup();
	// In the below, 'rng_type' and 'rng_r' are defined in FDTDvar.h      
	rng_type = gsl_rng_default;
	rng_r = gsl_rng_alloc (rng_type);
	gsl_rng_set(rng_r,1234);  // seed value = 1234 
	//////////////////////////////////////
	
	printf("parameters...ok\n");
}

void set_sigma_order(float oxl, float oxr, float oyl, float oyr, float ozl, float ozr)
{
	orderxl = oxl; orderyl = oyl; orderzl = ozl;  //default value = 3.5
	orderxr = oxr; orderyr = oyr; orderzr = ozr;  //default value = 3.5
}

void set_sigma_max(float axl, float axr, float ayl, float ayr, float azl, float azr)
{
	sig_axl = axl; sig_ayl = ayl; sig_azl = azl;   //default value = 1.0
	sig_axr = axr; sig_ayr = ayr; sig_azr = azr;   //default value = 1.0
}

void set_kappa(float kappa_x, float kappa_y, float kappa_z)
{
	kx = kappa_x;  //defaut value = 1.0
	ky = kappa_y;  //defaut value = 1.0
	kz = kappa_z;  //defaut value = 1.0
}

void Hz_parity(int x,int y,int z)
{
	xparity=x;
	yparity=y;
	zparity=z;

	if(xparity==1 || xparity==-1)
	{
		misize=cisize+2;
		pisize=cisize+1;
	}

	if(yparity==1 || yparity==-1)
	{
		mjsize=cjsize+2;
		pjsize=cjsize+1;
	}

	if(zparity==1 || zparity==-1)
	{
		mksize=cksize+2;
		pksize=cksize+1;
	}
}

void periodic_boundary(int x_on, int y_on, float k_x, float k_y)
{
	if( x_on == 1)
	{
		wave_vector_x = k_x;
		use_periodic_x = x_on; // 1:on, 0:off  default=0
		misize = isize+2; 
		pisize = isize+1;
		pmlil=-5;  //See the range of position[i][j][k] in memory.c
		pmlir=-5;
	}
	if( y_on == 1)
	{		
		wave_vector_y = k_y;
		use_periodic_y = y_on;
		mjsize = jsize+2; 
		pjsize = jsize+1;
		pmljl=-5;
		pmljr=-5;
	}
}

int lcm3(int a, int b, int c)
{
	return( lcm2(a, lcm2(b,c)) );
}

int lcm2(int m, int n)
{
	return( m*n/gcd2(m,n) );  //use the relation ; gcd(m,n)*lcm(m,n)=m*n
}

int gcd2(int m, int n)
{
	int temp;
	int r; //remainder

	if(n>m) //Make m>n
	{
		temp = m;
		m = n;
		n = temp;
	}
	while( (r= m%n) != 0)  //Euclid algorithm
	{
		m=n;
		n=r;
	}

	return(n);
}

float *float_1d_memory_for_parameter(int imax)
{
	int i;
	float *memory;

	memory=(float *)malloc(sizeof(float)*imax);
	for(i=0;i<imax;i++) memory[i]=0.0;

	return memory;
}