example.c

program FDTD for the photonic crystal structure

#include "~/bin/pFDTD.h"
#include "~/bin/pFDTDvar.h"

main()
/* This function is a main part of 3D FDTD code. */
{
	float x,y;
	int i,j,k;
	float bottom;
	float N; 
	float WC, *OMEGA, DETECT;
	float NROW, NA, Nfree;
	float R, T, Rm; 
	long DT,DD,DS;

//----------------------------------------------------------------------------//

	/////////////////////////////////////////
	/////// Basic Parameters (START) ////////
	/////////////////////////////////////////

	N = 12;  // in-plane domain size
	
	//bottom = 2.25;
	//shift = (4+bottom)/2 - 1.75;  
	shift = 0;

	///// dipole source parameter /////
	DT=500;
	DS=DT*6;
	DD=15000;
	WC=0.33;

	//----( far-field calculation )----//
	SpecN = 2;
	OMEGA = (float *)malloc(sizeof(float)*SpecN);
	///// set frequencies ////
	//for(i=0; i<SpecN; i++)
		//OMEGA[i] = 0.318+0.004*i;
	OMEGA[0] = 0.3179;
	OMEGA[1] = 0.3442;
	//////////////////////////
	DETECT = 0.27; 
	NROW = 512;
	NA = 0.85;
	Nfree = 1.0;
	//----------------------------------//

	//structure_size(N,N,4+bottom); 
	structure_size(N,N,3); 
	lattice_size(10,10,10);
	//non_uniform_grid("z",-0.2,0,30);
	pml_size(10,10,10,10,10,10);
	set_default_parameter(2);
	Hz_parity(1,1,-1);  // parities for Hz-field //
	memory();

	///////////////////////////////////////
	/////// Basic Parameters (END) ////////
	///////////////////////////////////////
//----------------------------------------------------------------------------//

	/////////////////////////////////////////
	/////// Input Structure (START) /////////
	/////////////////////////////////////////

	////// slab structure //////////
	R = 0.35;
	T = 0.5;
	Rm = 0.25;
	////////////////////////////////

	background(1.0);
	// dielectric slab
	input_object("block",EMP,0,0,0+shift,N,N,T,11.56);
	//input_object("block",EMP,0,0,1-(4+bottom)/2,N,N,2,11.56);

	// periodic lattice
	for(x=-(N/2);x<=(N/2);x++)
	for(y=-(N/2);y<=(N/2);y++)
	{
		input_object("rod",EMP,x,sqrt(3)*y,0+shift,R,T,0,1);
		input_object("rod",EMP,x-0.5,sqrt(3)*(y+0.5),0+shift,R,T,0,1);
	}

	// fill
	input_object("rod",EMP,0,0,0+shift,(R+0.01),T,0,11.56); //center
	input_object("rod",EMP,-1,0,0+shift,(R+0.01),T,0,11.56); //1
	input_object("rod",EMP,1,0,0+shift,(R+0.01),T,0,11.56); //2
	input_object("rod",EMP,0.5,sqrt(3)*0.5,0+shift,(R+0.01),T,0,11.56); //3
	input_object("rod",EMP,-0.5,sqrt(3)*0.5,0+shift,(R+0.01),T,0,11.56); //4
	input_object("rod",EMP,0.5,-sqrt(3)*0.5,0+shift,(R+0.01),T,0,11.56); //5
	input_object("rod",EMP,-0.5,-sqrt(3)*0.5,0+shift,(R+0.01),T,0,11.56); //6

	// dig
	input_object("rod",EMP,(1+R-Rm)*-1,0,0+shift,Rm,T,0,1); //1
	input_object("rod",EMP,(1+R-Rm)*1,0,0+shift,Rm,T,0,1); //2
	input_object("rod",EMP,(1+R-Rm)*0.5,1.1*sqrt(3)*0.5,0+shift,Rm,T,0,1); //3
	input_object("rod",EMP,(1+R-Rm)*-0.5,1.1*sqrt(3)*0.5,0+shift,Rm,T,0,1); //4
	input_object("rod",EMP,(1+R-Rm)*0.5,1.1*-sqrt(3)*0.5,0+shift,Rm,T,0,1); //5
	input_object("rod",EMP,(1+R-Rm)*-0.5,1.1*-sqrt(3)*0.5,0+shift,Rm,T,0,1); //6

	/*
	// real structure
	input_object("contour","mono_new",231,189,0,60,0.551,(20/30.615),1.0);
	input_object("contour","post_matrix",224,203,-1.25,10,2.5,(20/26.6)*(32.37/285)*0.55,11.56);
	*/

	// side blocks
	input_object("block",EMP,(N/2),0,0+shift,1.2,N+1,T,11.56);
	input_object("block",EMP,-(N/2),0,0+shift,1.2,N+1,T,11.56);
	input_object("block",EMP,0,(N/2),0+shift,N+1,1.2,T,11.56);
	input_object("block",EMP,0,-(N/2),0+shift,N+1,1.2,T,11.56);

	///////////////////////////////////////
	/////// Input Structure (END) /////////
	///////////////////////////////////////
//----------------------------------------------------------------------------//

	make_epsilon();
	out_epsilon("x",0,"epsilon.x");
	out_epsilon("y",0,"epsilon.y");
	out_epsilon("z",0+shift,"epsilon.z");

	coefficient(); 
//----------------------------------------------------------------------------//
	/////////////////////////////////////////////////////////////
	//////// Do something during the iterations (START) /////////
	/////////////////////////////////////////////////////////////	
	for(t=0;t<DD;t++)                        
	{	
		Gaussian_dipole_source("Hz",-0,-0.5,0+shift,WC,0,3*DT,DT);
		Gaussian_dipole_source("Hz",-0.2,-0.2,0+shift,WC,0,3*DT,DT);
		Gaussian_dipole_source("Hz",-0.4,-0.3,0+shift,WC,0,3*DT,DT);

		propagate();

		out_point("Hz",-0,-0.5,0+shift,0,DS,"source.dat");
		out_point("Hz",-0,-0.5,0+shift,DS,DD,"mode1.dat");

		if((DS+100)<t && t<DD) // farfield parameters
			far_field_param(OMEGA, DETECT+shift);

		/*
		if(DD-300<t && t<DD) // Qv and Qh calculation
		{	
			total_E_energy();
			total_E2();
			Poynting_total();
			Poynting_side(0.75,0);  //half width of the strip
		}

		if(DD-10<t && t<DD && t%2==0)
		{
			out_plane("Hz","z",0+shift,".Hz");
		}*/
	}

	///////////////////////////////////////////////////////////
	//////// Do something during the iterations (END) /////////
	///////////////////////////////////////////////////////////
//----------------------------------------------------------------------------//

	print_real_and_imag(0);
	print_real_and_imag_2n_size(NROW,0);

	far_field_FFT(NROW, NA, Nfree, OMEGA, 0);

	transform_farfield(NROW, 101, "rad_tot_", 0);
	transform_farfield(NROW, 101, "rad_Ex_", 0);
	transform_farfield(NROW, 101, "rad_Ey_", 0);
	transform_farfield(NROW, 101, "rad_Et_", 0);
	transform_farfield(NROW, 101, "rad_Ep_", 0);

	printf("Calculation Complete!\n");
}