farfield.c

program FDTD for the photonic crystal structure

	int m;
   	float *real,*imag;

   	/* Transform the rows */
   	real = (float *)malloc(nx * sizeof(float));
   	imag = (float *)malloc(nx * sizeof(float));

	m = pow_finder(nx);
   	for (j=0;j<ny;j++)
	{
      		for (i=0;i<nx;i++)
			{
         		real[i] = c[i][j].real;
         		imag[i] = c[i][j].imag;
      		}
			printf("i=%d,j=%d \r",i,j);
      		FFT1D(dir,m,real,imag);
      		for (i=0;i<nx;i++)
			{
         		c[i][j].real = real[i];
         		c[i][j].imag = imag[i];
      		}
   	}
	printf("\n");
   	free(real);
   	free(imag);

   	/* Transform the columns */
   	real = (float *)malloc(ny * sizeof(float));
   	imag = (float *)malloc(ny * sizeof(float));
   	
	m = pow_finder(ny);
   	for (i=0;i<nx;i++) 
	{
      		for (j=0;j<ny;j++)
			{
         		real[j] = c[i][j].real;
         		imag[j] = c[i][j].imag;
      		}
			printf("i=%d,j=%d \r",i,j);
      		FFT1D(dir,m,real,imag);
      		for (j=0;j<ny;j++) 
			{
         		c[i][j].real = real[j];
         		c[i][j].imag = imag[j];
      		}
   	}
	printf("\n");
   	free(real);
   	free(imag);

   	return(TRUE);
}

int FFT1D(int dir,int m,float *x,float *y)
{
   	long nn,i,i1,j,k,i2,l,l1,l2;
   	float c1,c2,tx,ty,t1,t2,u1,u2,z;

   	/* Calculate the number of points */
   	nn = 1;
   	for (i=0;i<m;i++)
      		nn *= 2;

   	/* Do the bit reversal */
   	i2 = nn >> 1;
   	j = 0;
   	for (i=0;i<nn-1;i++)
	{
      		if (i < j)
			{
         		tx = x[i];
         		ty = y[i];
         		x[i] = x[j];
         		y[i] = y[j];
         		x[j] = tx;
         		y[j] = ty;
      		}
      		k = i2;
      		while (k <= j)
			{
         		j -= k;
         		k >>= 1;
     	 	}
      		j += k;
   	}

   	/* Compute the FFT */
   	c1 = -1.0;
   	c2 = 0.0;
   	l2 = 1;
   	for (l=0;l<m;l++)
	{
      		l1 = l2;
      		l2 <<= 1;
      		u1 = 1.0;
      		u2 = 0.0;
      		for (j=0;j<l1;j++)
			{
         		for (i=j;i<nn;i+=l2)
				{
            			i1 = i + l1;
            			t1 = u1 * x[i1] - u2 * y[i1];
            			t2 = u1 * y[i1] + u2 * x[i1];
            			x[i1] = x[i] - t1;
            			y[i1] = y[i] - t2;
            			x[i] += t1;
           			y[i] += t2;
         		}
         		z =  u1 * c1 - u2 * c2;
         		u2 = u1 * c2 + u2 * c1;
         		u1 = z;
      		}
      		c2 = sqrt((1.0 - c1) / 2.0);
      		if (dir == 1)
         	c2 = -c2;
      		c1 = sqrt((1.0 + c1) / 2.0);
   	}

   	/* Scaling for forward transform */
   	if (dir == 1)
	{
      		for (i=0;i<nn;i++)
			{
         		x[i] /= (float)nn;
         		y[i] /= (float)nn;
      		}
   	}

   	return(TRUE);
}

int pow_finder(int N)
{
	int power=0;
	
	while( N != 1 )
	{
		N = (int)(N/2);
		power ++;
	}

	return(power);
}

void make_file_name(int mm)
{
	sprintf(name_freq,".ri%02d",mm);

	sprintf(name_Ex_real,"Ex_real");
	sprintf(name_Ex_imag,"Ex_imag");
	strcat(name_Ex_real,name_freq);
	strcat(name_Ex_imag,name_freq);

	sprintf(name_Ey_real,"Ey_real");
	sprintf(name_Ey_imag,"Ey_imag");
	strcat(name_Ey_real,name_freq);
	strcat(name_Ey_imag,name_freq);

	sprintf(name_Hx_real,"Hx_real");
	sprintf(name_Hx_imag,"Hx_imag");
	strcat(name_Hx_real,name_freq);
	strcat(name_Hx_imag,name_freq);

	sprintf(name_Hy_real,"Hy_real");
	sprintf(name_Hy_imag,"Hy_imag");
	strcat(name_Hy_real,name_freq);
	strcat(name_Hy_imag,name_freq);
	
	sprintf(FFT_E_int2,"FFT_E_int2_");
	sprintf(FFT_E_log,"FFT_E_log_");
	strcat(FFT_E_int2,name_freq);
	strcat(FFT_E_log,name_freq);
	sprintf(FFT_H_int2,"FFT_H_int2_");
	sprintf(FFT_H_log,"FFT_H_log_");
	strcat(FFT_H_int2,name_freq);
	strcat(FFT_H_log,name_freq);

	sprintf(radiation_tot,"rad_tot_");
	sprintf(radiation_Et,"rad_Et_");
	sprintf(radiation_Ep,"rad_Ep_");
	sprintf(radiation_Ex,"rad_Ex_");
	sprintf(radiation_Ey,"rad_Ey_");
	strcat(radiation_tot,name_freq);
	strcat(radiation_Et,name_freq);
	strcat(radiation_Ep,name_freq);
	strcat(radiation_Ex,name_freq);
	strcat(radiation_Ey,name_freq);

	sprintf(polar,"polarizaton_");
	sprintf(P_tot_name,"P_tot_");
	sprintf(P_the_name,"P_the_");
	sprintf(P_phi_name,"P_phi_");
	strcat(polar,name_freq);
	strcat(P_tot_name,name_freq);
	strcat(P_the_name,name_freq);
	strcat(P_phi_name,name_freq);
}

void set_global_variable(int NROW, float OMEGA, float Nfree, int *POLROT, int *dir, float *Eta, float *k)
{
	*POLROT = 36;
	*dir = 1; // forward transformation
	*Eta = 377/Nfree;
	*k = floor(NROW*OMEGA/lattice_x);
	printf("radius k=%4.1f\n",*k);

	printf("Now... Caculating far-field pattern......\n");
}

//////////////////////////////
/// Reading source         ///
//////////////////////////////

void reading_source_Ex(int row, int col, float **efieldx_real, float **efieldx_imag)
{
	int i, j; 
	FILE *stream;

	stream = fopen(name_Ex_real,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				efieldx_real[i][j] = 0.0;
			else
				efieldx_real[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Ex_imag,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				efieldx_imag[i][j] = 0.0;
			else
				efieldx_imag[i][j] = atof(string);
		}
	}
	fclose(stream);

	printf("reading Ex check...\n");
}

void reading_source_Ey(int row, int col, float **efieldy_real, float **efieldy_imag)
{
	int i, j; 
	FILE *stream;

	stream = fopen(name_Ey_real,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				efieldy_real[i][j] = 0.0;
			else
				efieldy_real[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Ey_imag,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				efieldy_imag[i][j] = 0.0;
			else
				efieldy_imag[i][j] = atof(string);
		}
	}
	fclose(stream);

	printf("reading Ey check...\n");
}

void reading_source_Hx(int row, int col, float **hfieldx_real, float **hfieldx_imag)
{
	int i, j; 
	FILE *stream;

	stream = fopen(name_Hx_real,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				hfieldx_real[i][j] = 0.0;
			else
				hfieldx_real[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Hx_imag,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				hfieldx_imag[i][j] = 0.0;
			else
				hfieldx_imag[i][j] = atof(string);
		}
	}
	fclose(stream);

	printf("reading Hx check...\n");
}

void reading_source_Hy(int row, int col, float **hfieldy_real, float **hfieldy_imag)
{
	int i, j; 
	FILE *stream;

	stream = fopen(name_Hy_real,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				hfieldy_real[i][j] = 0.0;
			else
				hfieldy_real[i][j] = atof(string);
		}
	}
	fclose(stream);

	stream = fopen(name_Hy_imag,"rt");
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			fscanf(stream, "%s", string);
			if( strcmp("nan",string) == 0 )
				hfieldy_imag[i][j] = 0.0;
			else
				hfieldy_imag[i][j] = atof(string);
		}
	}
	fclose(stream);

	printf("reading Hy check...\n");
}

//////////////////////////////
/// Fourier Transformation ///
//////////////////////////////

void fourier_transform_Nx_Hy(int row, int col, int dir, float **Nx_real, float **Nx_imag, float **hfieldy_real, float **hfieldy_imag)
{
	int i, j; 
	struct COMPLEX **FFT; // FFT input and output
	int FFT_return;

	FFT = (struct COMPLEX **)malloc(sizeof(struct COMPLEX *)*col);
	for(i=0; i<col; i++)
		FFT[i] = (struct COMPLEX *)malloc(sizeof(struct COMPLEX)*row);
	printf("FFT allocation complete....\n");

	///------------( Nx calculation )------------
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			FFT[i][j].real = hfieldy_real[i][j];
			FFT[i][j].imag = hfieldy_imag[i][j];
		}
	}	
	printf("starting Nx caculation (FFT)...\n");
	FFT_return = FFT2D(FFT, col, row, dir);
	printf("FFT check... (FFT_return=%d)\n",FFT_return);

	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			Nx_real[i][j] = -FFT[i][j].real;
			Nx_imag[i][j] = -FFT[i][j].imag;
		}
	}

	for(i=0; i<col; i++)
	{
		free(FFT[i]);
	}
	free(FFT);
}

void fourier_transform_Ny_Hx(int row, int col, int dir, float **Ny_real, float **Ny_imag, float **hfieldx_real, float **hfieldx_imag)
{
	int i, j; 
	struct COMPLEX **FFT; // FFT input and output
	int FFT_return;

	FFT = (struct COMPLEX **)malloc(sizeof(struct COMPLEX *)*col);
	for(i=0; i<col; i++)
		FFT[i] = (struct COMPLEX *)malloc(sizeof(struct COMPLEX)*row);
	printf("FFT allocation complete....\n");

	///------------( Ny calculation )------------
	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			FFT[i][j].real = hfieldx_real[i][j];
			FFT[i][j].imag = hfieldx_imag[i][j];
		}
	}
	printf("starting Ny caculation (FFT)...\n");
	FFT_return = FFT2D(FFT, col, row, dir);
	printf("FFT check... (FFT_return=%d)\n",FFT_return);

	for(j=row-1; j>=0; j--)
	{
		for(i=0; i<col; i++)
		{
			Ny_real[i][j] = FFT[i][j].real;
			Ny_imag[i][j] = FFT[i][j].imag;
		}
	}