farfield.c

program FDTD for the photonic crystal structure

	half_col = (int)(col/2);
	half_row = (int)(row/2);

	row_s = (int)(10*Nfree*k);
	col_s = (int)(10*Nfree*k);

	half_col_s = (int)(col_s/2);
	half_row_s = (int)(row_s/2);

	for(j=(half_row)+(half_row_s); j>=(half_row)-(half_row_s); j--)
	{
		for(i=(half_col)-(half_col_s); i<=(half_col)+(half_col_s); i++)
		{
			Ly_s_real[i-(half_col)+(half_col_s)][j-(half_row)+(half_row_s)] = Ly_real[i][j];
			Ly_s_imag[i-(half_col)+(half_col_s)][j-(half_row)+(half_row_s)] = Ly_imag[i][j];
		}
	}
}

//////////////////////////////
/// FFT E & FFT H          ///
//////////////////////////////

void FT_data_writing(float **Nx_s_real, float **Nx_s_imag, float **Ny_s_real, float **Ny_s_imag, float **Lx_s_real, float **Lx_s_imag, float **Ly_s_real, float **Ly_s_imag, int row_s, int col_s, float Nfree, float k)
{
	int half_row_s, half_col_s;
	float FFT_Emax, FFT_Hmax;  // FFT normalization
	float temp;
	FILE *stream;
	int i,j;

	half_col_s = (int)(col_s/2);
	half_row_s = (int)(row_s/2);

	FFT_Emax = 0.0;
	FFT_Hmax = 0.0;

	//////////[ FFT(Ex)^2 + FFT(Ey)^2 ]///////////
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<=(col_s); i++)
		{
				temp = Lx_s_real[i][j]*Lx_s_real[i][j]+Ly_s_real[i][j]*Ly_s_real[i][j]+Lx_s_imag[i][j]*Lx_s_imag[i][j]+Ly_s_imag[i][j]*Ly_s_imag[i][j];
				if( temp > FFT_Emax)
					FFT_Emax = temp;
		}
	}

	stream = fopen(FFT_E_int2, "wt");
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<=(col_s); i++)
		{
			fprintf(stream,"%g\t",(Lx_s_real[i][j]*Lx_s_real[i][j]+Ly_s_real[i][j]*Ly_s_real[i][j]+Lx_s_imag[i][j]*Lx_s_imag[i][j]+Ly_s_imag[i][j]*Ly_s_imag[i][j])/FFT_Emax);
		}
		fprintf(stream,"\n");
	}
	fclose(stream);

	stream = fopen(FFT_E_log, "wt");
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<=(col_s); i++)
		{
			fprintf(stream,"%g\t",log10((Lx_s_real[i][j]*Lx_s_real[i][j]+Ly_s_real[i][j]*Ly_s_real[i][j]+Lx_s_imag[i][j]*Lx_s_imag[i][j]+Ly_s_imag[i][j]*Ly_s_imag[i][j])/FFT_Emax));
		}
		fprintf(stream,"\n");
	}
	fclose(stream);

	//////////[ FFT(Hx)^2 + FFT(Hy)^2 ]///////////
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<=(col_s); i++)
		{
				temp = Nx_s_real[i][j]*Nx_s_real[i][j]+Ny_s_real[i][j]*Ny_s_real[i][j]+Nx_s_imag[i][j]*Nx_s_imag[i][j]+Ny_s_imag[i][j]*Ny_s_imag[i][j];
				if( temp > FFT_Hmax)
					FFT_Hmax = temp;
		}
	}

	stream = fopen(FFT_H_int2, "wt");
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<=(col_s); i++)
		{
			fprintf(stream,"%g\t",(Nx_s_real[i][j]*Nx_s_real[i][j]+Ny_s_real[i][j]*Ny_s_real[i][j]+Nx_s_imag[i][j]*Nx_s_imag[i][j]+Ny_s_imag[i][j]*Ny_s_imag[i][j])/FFT_Hmax);
		}
		fprintf(stream,"\n");
	}
	fclose(stream);

	stream = fopen(FFT_H_log, "wt");
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<=(col_s); i++)
		{
			fprintf(stream,"%g\t",log10((Nx_s_real[i][j]*Nx_s_real[i][j]+Ny_s_real[i][j]*Ny_s_real[i][j]+Nx_s_imag[i][j]*Nx_s_imag[i][j]+Ny_s_imag[i][j]*Ny_s_imag[i][j])/FFT_Hmax));
		}
		fprintf(stream,"\n");
	}
	fclose(stream);
}	

void set_radiation_variable(float Nfree, float k, float Eta, float **Nx_s_real, float **Nx_s_imag, float **Ny_s_real, float **Ny_s_imag, float **Lx_s_real, float **Lx_s_imag, float **Ly_s_real, float **Ly_s_imag, float **Nt_s_real, float **Nt_s_imag, float **Np_s_real, float **Np_s_imag, float **Lt_s_real, float **Lt_s_imag, float **Lp_s_real, float **Lp_s_imag, float **cosp, float **sinp, float **cost, float **ex_real, float **ex_imag, float **ey_real, float **ey_imag)
{
	int i,j;
	int row_s, col_s;           //shrinked row & col
	int half_row_s, half_col_s;

	row_s = (int)(10*Nfree*k);
	col_s = (int)(10*Nfree*k);

	half_col_s = (int)(col_s/2);
	half_row_s = (int)(row_s/2);

	printf("calculating N(theta), N(phi).....\n");

	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k )
			{
				cosp[i][j] = (i-half_col_s)/sqrt( (i-half_col_s)*(i-half_col_s) + (j-half_row_s)*(j-half_row_s) );
				sinp[i][j] = (j-half_row_s)/sqrt( (i-half_col_s)*(i-half_col_s) + (j-half_row_s)*(j-half_row_s) );
				cost[i][j] = sqrt( 1 - ((i-half_col_s)*(i-half_col_s) + (j-half_row_s)*(j-half_row_s))/(Nfree*Nfree*k*k) );
			}
			else
			{
				cosp[i][j] = 0.0;
				sinp[i][j] = 0.0;
				cost[i][j] = 0.0;
			}
		}
	}
	cosp[half_col_s][half_row_s] = 1/sqrt(2);
	sinp[half_col_s][half_row_s] = 1/sqrt(2);
	cost[half_col_s][half_row_s] = 1;

	printf("trigonometric fucntion define...\n");

	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			Nt_s_real[i][j] = (Nx_s_real[i][j]*cosp[i][j]+Ny_s_real[i][j]*sinp[i][j])*cost[i][j];
			Nt_s_imag[i][j] = (Nx_s_imag[i][j]*cosp[i][j]+Ny_s_imag[i][j]*sinp[i][j])*cost[i][j];
			Np_s_real[i][j] = (-Nx_s_real[i][j]*sinp[i][j]+Ny_s_real[i][j]*cosp[i][j]);
			Np_s_imag[i][j] = (-Nx_s_imag[i][j]*sinp[i][j]+Ny_s_imag[i][j]*cosp[i][j]);
			Lt_s_real[i][j] = (Lx_s_real[i][j]*cosp[i][j]+Ly_s_real[i][j]*sinp[i][j])*cost[i][j];
			Lt_s_imag[i][j] = (Lx_s_imag[i][j]*cosp[i][j]+Ly_s_imag[i][j]*sinp[i][j])*cost[i][j];
			Lp_s_real[i][j] = (-Lx_s_real[i][j]*sinp[i][j]+Ly_s_real[i][j]*cosp[i][j]);
			Lp_s_imag[i][j] = (-Lx_s_imag[i][j]*sinp[i][j]+Ly_s_imag[i][j]*cosp[i][j]);
		}
	}
	printf("Nt, Np, Lt, Lp ok!.....\n");

	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			ex_real[i][j] = -(Np_s_real[i][j] - Lt_s_real[i][j]/Eta)*sinp[i][j] + (Nt_s_real[i][j] + Lp_s_real[i][j]/377)*cost[i][j]*cosp[i][j];
			ex_imag[i][j] = -(Np_s_imag[i][j] - Lt_s_imag[i][j]/Eta)*sinp[i][j] + (Nt_s_imag[i][j] + Lp_s_imag[i][j]/377)*cost[i][j]*cosp[i][j];
			ey_real[i][j] = (Np_s_real[i][j] - Lt_s_real[i][j]/Eta)*cosp[i][j] + (Nt_s_real[i][j] + Lp_s_real[i][j]/Eta)*cost[i][j]*sinp[i][j];
			ey_imag[i][j] = (Np_s_imag[i][j] - Lt_s_imag[i][j]/Eta)*cosp[i][j] + (Nt_s_imag[i][j] + Lp_s_imag[i][j]/Eta)*cost[i][j]*sinp[i][j];
		}
	}
	printf("Ex, Ey ok!......\n");
}

void calc_radiation(float Nfree, float k, float Eta, float **Nt_s_real, float **Nt_s_imag, float **Np_s_real, float **Np_s_imag, float **Lt_s_real, float **Lt_s_imag, float **Lp_s_real, float **Lp_s_imag, float **cosp, float **sinp, float **cost, float **ex_real, float **ex_imag, float **ey_real, float **ey_imag, float **radint, float **etheta, float **ephi, float **eintx, float **einty, float *P_tot, float *P_the, float *P_phi)
{
	int i,j;
	int row_s, col_s;           //shrinked row & col
	int half_row_s, half_col_s;

	row_s = (int)(10*Nfree*k);
	col_s = (int)(10*Nfree*k);

	half_col_s = (int)(col_s/2);
	half_row_s = (int)(row_s/2);

	///------------( total intensity )------------
	printf("calculating radiation intensity.....\n");

	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				radint[i][j] = ((Nt_s_real[i][j]+Lp_s_real[i][j]/Eta)*(Nt_s_real[i][j]+Lp_s_real[i][j]/Eta)+ (Nt_s_imag[i][j]+Lp_s_imag[i][j]/Eta)*(Nt_s_imag[i][j]+Lp_s_imag[i][j]/Eta)+(Np_s_real[i][j]-Lt_s_real[i][j]/Eta)*(Np_s_real[i][j]-Lt_s_real[i][j]/Eta)+(Np_s_imag[i][j]-Lt_s_imag[i][j]/Eta)*(Np_s_imag[i][j]-Lt_s_imag[i][j]/Eta));
			else
				radint[i][j] = 0.0;
		}
	}
	printf("radint ok!....\n");

	///------------( Total Radiated Power )------------
	*P_tot = 0.0; //initialization
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				*P_tot = *P_tot + radint[i][j]/cost[i][j];
		}
	}


	///------------( Etheta intensity )------------
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				etheta[i][j] = ((Nt_s_real[i][j]+Lp_s_real[i][j]/Eta)*(Nt_s_real[i][j]+Lp_s_real[i][j]/Eta)
					+ (Nt_s_imag[i][j]+Lp_s_imag[i][j]/Eta)*(Nt_s_imag[i][j]+Lp_s_imag[i][j]/Eta));
			else
				etheta[i][j] = 0;
		}
	}

	///------------( Etheta Radiated Power )------------
	*P_the = 0.0; //initialization
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				*P_the = *P_the + etheta[i][j]/cost[i][j];
		}
	}

	///------------( Ephi intensity )------------
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				ephi[i][j] = ((Np_s_real[i][j]-Lt_s_real[i][j]/Eta)*(Np_s_real[i][j]-Lt_s_real[i][j]/Eta)
					+ (Np_s_imag[i][j]-Lt_s_imag[i][j]/Eta)*(Np_s_imag[i][j]-Lt_s_imag[i][j]/Eta));
			else
				ephi[i][j] = 0;
		}
	}

	///------------( Ephi Radiated Power )------------
	*P_phi = 0.0; //initialization
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				*P_phi = *P_phi + ephi[i][j]/cost[i][j];
		}
	}
	printf("etheta, ephi  ok!...\n");

	///------------( Ex intensity )------------
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				eintx[i][j] = (ex_real[i][j]*ex_real[i][j] + ex_imag[i][j]*ex_imag[i][j]);
			else
				eintx[i][j] = 0;
		}
	}

	///------------( Ey intensity )------------
	for(j=row_s; j>=0; j--)
	{
		for(i=0; i<(col_s+1); i++)
		{
			if( (i-half_col_s)*(i-half_col_s)+(j-half_row_s)*(j-half_row_s) < Nfree*Nfree*k*k) 
				einty[i][j] = (ey_real[i][j]*ey_real[i][j] + ey_imag[i][j]*ey_imag[i][j]);
			else
				einty[i][j] = 0;
		}
	}
	printf("eintx, einty ok!...\n");
}

void print_radiation(float Nfree, float k, float Eta, float **radint, float **etheta, float **ephi, float **eintx, float **einty, float *P_tot, float *P_the, float *P_phi)
{
	int i,j;
	int row_s, col_s;           //shrinked row & col
	int half_row_s, half_col_s;
	FILE *stream;

	row_s = (int)(10*Nfree*k);
	col_s = (int)(10*Nfree*k);

	half_col_s = (int)(col_s/2);
	half_row_s = (int)(row_s/2);

	///------------( total intensity )------------
	printf("radiation result writing......\n");
	stream = fopen(radiation_tot, "wt");
	for(j=half_row_s+(int)(Nfree*k)-1; j>=half_row_s-(int)(Nfree*k)+1; j--)
	{
		for(i=half_col_s-(int)(Nfree*k)+1; i<=half_col_s+(int)(Nfree*k)-1; i++)
		{
				fprintf(stream,"%g\t",radint[i][j]);
		}
		fprintf(stream,"\n");
	}
	fclose(stream);

	///------------( Total Radiated Power )------------
	printf("radiated power writing.......\n");
	stream = fopen(P_tot_name,"wt");
	fprintf(stream,"%g",P_tot);
	fclose(stream);

	///------------( Etheta Radiated Power )------------
	printf("radiated power writing.......\n");
	stream = fopen(P_the_name,"wt");
	fprintf(stream,"%g",P_the);
	fclose(stream);

	///------------( Ephi Radiated Power )------------
	printf("radiated power writing.......\n");
	stream = fopen(P_phi_name,"wt");
	fprintf(stream,"%g",P_phi);
	fclose(stream);

	///------------( Etheta intensity )------------
	stream = fopen(radiation_Et, "wt");
	for(j=half_row_s+(int)(Nfree*k)-1; j>=half_row_s-(int)(Nfree*k)+1; j--)
	{
		for(i=half_col_s-(int)(Nfree*k)+1; i<=half_col_s+(int)(Nfree*k)-1; i++)
		{
			fprintf(stream,"%g\t",etheta[i][j]);
		}
		fprintf(stream,"\n");
	}
	fclose(stream);

	///------------( Ephi intensity )------------
	stream = fopen(radiation_Ep, "wt");
	for(j=half_row_s+(int)(Nfree*k)-1; j>=half_row_s-(int)(Nfree*k)+1; j--)
	{
		for(i=half_col_s-(int)(Nfree*k)+1; i<=half_col_s+(int)(Nfree*k)-1; i++)
		{
			fprintf(stream,"%g\t",ephi[i][j]);
		}
		fprintf(stream,"\n");
	}
	fclose(stream);

	///------------( Ex intensity )------------
	stream = fopen(radiation_Ex, "wt");
	for(j=half_row_s+(int)(Nfree*k)-1; j>=half_row_s-(int)(Nfree*k)+1; j--)
	{
		for(i=half_col_s-(int)(Nfree*k)+1; i<=half_col_s+(int)(Nfree*k)-1; i++)
		{
			fprintf(stream,"%g\t",eintx[i][j]);
		}
		fprintf(stream,"\n");
	}
	fclose(stream);

	///------------( Ey intensity )------------
	stream = fopen(radiation_Ey, "wt");
	for(j=half_row_s+(int)(Nfree*k)-1; j>=half_row_s-(int)(Nfree*k)+1; j--)
	{
		for(i=half_col_s-(int)(Nfree*k)+1; i<=half_col_s+(int)(Nfree*k)-1; i++)
		{
			fprintf(stream,"%g\t",einty[i][j]);
		}
		fprintf(stream,"\n");
	}
	fclose(stream);
}

/*
void calc_polar(float NA, float Nfree)
{
	Eu_real = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Eu_real[i] = (float *)malloc(sizeof(float)*(row+1));
	Eu_imag = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Eu_imag[i] = (float *)malloc(sizeof(float)*(row+1));
	Eu_int = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Eu_int[i] = (float *)malloc(sizeof(float)*(row+1));
	polarization = (float *)malloc(sizeof(float)*POLROT);

	for(l=0; l<POLROT; l++)
	{
		sum_Eu = 0.0;
		////// projection into u(phi) //////
		for(j=row; j>=0; j--)
		{
			for(i=0; i<(col+1); i++)
			{
				if( (i-half_col)*(i-half_col)+(j-half_row)*(j-half_row) < NA*NA*Nfree*Nfree*k*k)
				{
					Eu_real[i][j] = ex_real[i][j]*cos(2*pi*l/POLROT) + ey_real[i][j]*sin(2*pi*l/POLROT);
					Eu_imag[i][j] = ex_imag[i][j]*cos(2*pi*l/POLROT) + ey_imag[i][j]*sin(2*pi*l/POLROT);
					Eu_int[i][j] = Eu_real[i][j]*Eu_real[i][j] + Eu_imag[i][j]*Eu_imag[i][j];
					sum_Eu = sum_Eu + Eu_int[i][j]*sqrt(1-cost[i][j]*cost[i][j])/FFT_Emax;
				}
			}
		}
		printf("polarization [%d] = %g\n", l,sum_Eu);
		polarization[l] = sum_Eu;
	}

	///// writing file ///////
	stream = fopen(polar,"wt");
	for(l=0; l<POLROT; l++)
	{
		fprintf(stream,"%d\t%g\t\n",l*10,polarization[l]);
	}
	fprintf(stream,"%d\t%g\t\n",l*10,polarization[0]);
	fclose(stream);

	free(ex_real); free(ey_real); free(ex_imag); free(ey_imag);
	free(cosp); free(sinp); free(cost);
	free(radint); free(etheta); free(ephi);
	free(eintx); free(einty);
	free(Eu_real); free(Eu_imag); free(Eu_int);
	free(polarization);
}
*/