farfield.c

program FDTD for the photonic crystal structure

#include "./pFDTD.h"

/////////////////////////////////////////////
///////// Define global variables ///////////
/////////////////////////////////////////////
static char string[80];  // global name

//static int l;

//static int m,n;
//static char ch;

//static float **Eu_real, **Eu_imag, **Eu_int;
//static float *polarization;
//static float sum_Eu;
//static float Ex_delta, Ey_delta, Hx_delta, Hy_delta;

static char name_freq[10];
static char name_Ex_real[20], name_Ex_imag[20], name_Ey_real[20], name_Ey_imag[20];
static char name_Hx_real[20], name_Hx_imag[20], name_Hy_real[20], name_Hy_imag[20];
static char FFT_E_int2[20], FFT_E_log[20], FFT_H_int2[20], FFT_H_log[20];
static char radiation_tot[20], radiation_Et[20], radiation_Ep[20], radiation_Ex[20], radiation_Ey[20];
static char polar[20], P_tot_name[20], P_the_name[20], P_phi_name[20];
//--------------------------------------------------------------------//

int FFT2D(struct COMPLEX **c,int nx,int ny,int dir);
int FFT1D(int dir,int m,float *x,float *y);
int pow_finder(int N);
void make_file_name(int mm);
void set_global_variable(int NROW, float OMEGA, float Nfree, int *POLROT, int *dir, float *Eta, float *k);

//-----------------------------------------------------------
// reading (NROW)x(NROW) source data 
// then, store them in efieldx_real[][], efieldx_imag[][] ... 
//-----------------------------------------------------------
void reading_source_Ex(int row, int col, float **efieldx_real, float **efieldx_imag);
void reading_source_Ey(int row, int col, float **efieldy_real, float **efieldy_imag);
void reading_source_Hx(int row, int col, float **hfieldx_real, float **hfieldx_imag);
void reading_source_Hy(int row, int col, float **hfieldy_real, float **hfieldy_imag);

//------------------------------------------------------------
// Fourier transforming Ex, Ey, Hx, Hy to get Lx, Ly, Nx, Ny
//------------------------------------------------------------
void fourier_transform_Ly_Ex(int row, int col, int dir, float **Ly_real, float **Ly_imag, float **efieldx_real, float **efieldx_imag);
void fourier_transform_Lx_Ey(int row, int col, int dir, float **Lx_real, float **Lx_imag, float **efieldy_real, float **efieldy_imag);
void fourier_transform_Ny_Hx(int row, int col, int dir, float **Ny_real, float **Ny_imag, float **hfieldx_real, float **hfieldx_imag);
void fourier_transform_Nx_Hy(int row, int col, int dir, float **Nx_real, float **Nx_imag, float **hfieldy_real, float **hfieldy_imag);

//------------------------------
// Data shifting Lx,Ly,Nx,Ny
//------------------------------
void data_shifting_Ly(int row, int col, float **Ly_real, float **Ly_imag);
void data_shifting_Lx(int row, int col, float **Lx_real, float **Lx_imag);
void data_shifting_Ny(int row, int col, float **Ny_real, float **Ny_imag);
void data_shifting_Nx(int row, int col, float **Nx_real, float **Nx_imag);

//------------------------------------------------------------------------
// Make small size data (-5*k ~ +5*k) ; k is the radius of the light-cone 
//------------------------------------------------------------------------ 
void Ly_data_shrink(int row, int col, float **Ly_real, float **Ly_imag, float **Ly_s_real, float **Ly_s_imag, float Nfree, float k);
void Lx_data_shrink(int row, int col, float **Lx_real, float **Lx_imag, float **Lx_s_real, float **Lx_s_imag, float Nfree, float k);
void Nx_data_shrink(int row, int col, float **Nx_real, float **Nx_imag, float **Nx_s_real, float **Nx_s_imag, float Nfree, float k);
void Ny_data_shrink(int row, int col, float **Ny_real, float **Ny_imag, float **Ny_s_real, float **Ny_s_imag, float Nfree, float k);

//----------------------------------------------------------
// Calculate Nt,Np,Lt,Lp, ex_real,ex_imag, ey_real,ey_imag
//----------------------------------------------------------
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);

//------------------------------------------------
// Calculate radint, etheta, ephi, eintx, einty 
//------------------------------------------------
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);

//------------------------------------------------
// Print radint, etheta, ephi, eintx, einty 
//------------------------------------------------
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);

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);
//void calc_polar(float NA, float Nfree);

void far_field_param(float *OMEGA, float DETECT)
{
	int i, j, k;
	int mm;

	k = non_uniform_z_to_i(DETECT);

	for(j=1; j<=jsize-2; j++)
	{
		for(i=1; i<=isize-2; i++)
		{
			for(mm=0; mm<SpecN; mm++)
			{
				Ex_cos[i][j][mm]+=grid_value("Ex",i,j,k)*cos(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
				Ex_sin[i][j][mm]+=grid_value("Ex",i,j,k)*sin(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
				Ey_cos[i][j][mm]+=grid_value("Ey",i,j,k)*cos(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
				Ey_sin[i][j][mm]+=grid_value("Ey",i,j,k)*sin(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
				Hx_cos[i][j][mm]+=grid_value("Hx",i,j,k)*cos(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
				Hx_sin[i][j][mm]+=grid_value("Hx",i,j,k)*sin(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
				Hy_cos[i][j][mm]+=grid_value("Hy",i,j,k)*cos(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
				Hy_sin[i][j][mm]+=grid_value("Hy",i,j,k)*sin(2*pi*OMEGA[mm]*t/S_factor/ds_x/lattice_x);
			}
		}
	}	
}

void far_field_FFT(int NROW, float NA, float Nfree, float *OMEGA, int mm)
{
	int i;
	int row, col; // row num and column num
	int row_s, col_s; // shrinked row & col num

	int POLROT;
	int dir; // direction of the FFT
	float Eta;
	float k; // normalized k value, radius of the light-cone

	//// source matrix ////
	float **hfieldx_real, **hfieldx_imag;
	float **hfieldy_real, **hfieldy_imag;
	float **efieldx_real, **efieldx_imag;
	float **efieldy_real, **efieldy_imag;

	//// after FFT ////
	float **Nx_real, **Nx_imag, **Ny_real, **Ny_imag;
	float **Lx_real, **Lx_imag, **Ly_real, **Ly_imag;
	float **Nx_s_real, **Nx_s_imag, **Ny_s_real, **Ny_s_imag;   // s : shrinked size
	float **Lx_s_real, **Lx_s_imag, **Ly_s_real, **Ly_s_imag;

	float **Nt_s_real, **Nt_s_imag, **Np_s_real, **Np_s_imag;   // s : shrinked size
	float **Lt_s_real, **Lt_s_imag, **Lp_s_real, **Lp_s_imag;

	float **ex_real, **ex_imag, **ey_real, **ey_imag;  // shrinked size 

	float **cosp, **sinp, **cost;  // shrinked size

	float **radint; // radiation intensity
	float **etheta, **ephi;
	float **eintx, **einty;
	float P_tot, P_the, P_phi; // integrating the radiated power

	//------------------------------------------------------
	col = NROW, row = NROW;

	make_file_name(mm);
	set_global_variable(NROW, OMEGA[mm], Nfree, &POLROT, &dir, &Eta, &k);

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

	//////////////////////
	//---- Ex field ----//
	//////////////////////
	efieldx_real = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		efieldx_real[i] = (float *)malloc(sizeof(float)*row);
	efieldx_imag = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		efieldx_imag[i] = (float *)malloc(sizeof(float)*row);
	//------------------------------------------------------------
	reading_source_Ex(row, col, efieldx_real, efieldx_imag);
	//------------------------------------------------------------
	Ly_real = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Ly_real[i] = (float *)malloc(sizeof(float)*(row+1));
	Ly_imag = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Ly_imag[i] = (float *)malloc(sizeof(float)*(row+1));
	//-----------------------------------------------------------------------------------
	fourier_transform_Ly_Ex(row, col, dir, Ly_real, Ly_imag, efieldx_real, efieldx_imag);
	//-----------------------------------------------------------------------------------
	for(i=0; i<col; i++)
	{
		free(efieldx_real[i]); free(efieldx_imag[i]);
	}
		free(efieldx_real); free(efieldx_imag);
	//-------------------------------------------
	data_shifting_Ly(row, col, Ly_real, Ly_imag);
	//-------------------------------------------
	Ly_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Ly_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Ly_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Ly_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	//-------------------------------------------------------------------------
	Ly_data_shrink(row, col, Ly_real, Ly_imag, Ly_s_real, Ly_s_imag, Nfree, k);
	//-------------------------------------------------------------------------
	for(i=0; i<col+1; i++)
	{
		free(Ly_real[i]); free(Ly_imag[i]);
	}
		free(Ly_real); free(Ly_imag);

	//////////////////////
	//---- Ey field ----//
	//////////////////////
	efieldy_real = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		efieldy_real[i] = (float *)malloc(sizeof(float)*row);
	efieldy_imag = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		efieldy_imag[i] = (float *)malloc(sizeof(float)*row);
	//------------------------------------------------------------
	reading_source_Ey(row, col, efieldy_real, efieldy_imag);
	//------------------------------------------------------------
	Lx_real = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Lx_real[i] = (float *)malloc(sizeof(float)*(row+1));
	Lx_imag = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Lx_imag[i] = (float *)malloc(sizeof(float)*(row+1));
	//-----------------------------------------------------------------------------------
	fourier_transform_Lx_Ey(row, col, dir, Lx_real, Lx_imag, efieldy_real, efieldy_imag);
	//-----------------------------------------------------------------------------------
	for(i=0; i<col; i++)
	{
		free(efieldy_real[i]); free(efieldy_imag[i]);
	}
		free(efieldy_real); free(efieldy_imag);
	//-------------------------------------------
	data_shifting_Lx(row, col, Lx_real, Lx_imag);
	//-------------------------------------------
	Lx_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Lx_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Lx_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Lx_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	//-------------------------------------------------------------------------
	Lx_data_shrink(row, col, Lx_real, Lx_imag, Lx_s_real, Lx_s_imag, Nfree, k);
	//-------------------------------------------------------------------------
	for(i=0; i<col+1; i++)
	{
		free(Lx_real[i]); free(Lx_imag[i]);
	}
		free(Lx_real); free(Lx_imag);

	//////////////////////
	//---- Hx field ----//
	//////////////////////
	hfieldx_real = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		hfieldx_real[i] = (float *)malloc(sizeof(float)*row);
	hfieldx_imag = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		hfieldx_imag[i] = (float *)malloc(sizeof(float)*row);
	//------------------------------------------------------------
	reading_source_Hx(row, col, hfieldx_real, hfieldx_imag);
	//------------------------------------------------------------
	Ny_real = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Ny_real[i] = (float *)malloc(sizeof(float)*(row+1));
	Ny_imag = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Ny_imag[i] = (float *)malloc(sizeof(float)*(row+1));
	//-----------------------------------------------------------------------------------
	fourier_transform_Ny_Hx(row, col, dir, Ny_real, Ny_imag, hfieldx_real, hfieldx_imag);
	//-----------------------------------------------------------------------------------
	for(i=0; i<col; i++)
	{
		free(hfieldx_real[i]); free(hfieldx_imag[i]);
	}
		free(hfieldx_real); free(hfieldx_imag);
	//-------------------------------------------
	data_shifting_Ny(row, col, Ny_real, Ny_imag);
	//-------------------------------------------
	Ny_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Ny_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Ny_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Ny_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	//-------------------------------------------------------------------------
	Ny_data_shrink(row, col, Ny_real, Ny_imag, Ny_s_real, Ny_s_imag, Nfree, k);
	//-------------------------------------------------------------------------
	for(i=0; i<col+1; i++)
	{
		free(Ny_real[i]); free(Ny_imag[i]);
	}
		free(Ny_real); free(Ny_imag);

	//////////////////////
	//---- Hy field ----//
	//////////////////////
	hfieldy_real = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		hfieldy_real[i] = (float *)malloc(sizeof(float)*row);
	hfieldy_imag = (float **)malloc(sizeof(float *)*col);
	for(i=0; i<col; i++)
		hfieldy_imag[i] = (float *)malloc(sizeof(float)*row);
	//------------------------------------------------------------
	reading_source_Hy(row, col, hfieldy_real, hfieldy_imag);
	//------------------------------------------------------------
	Nx_real = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Nx_real[i] = (float *)malloc(sizeof(float)*(row+1));
	Nx_imag = (float **)malloc(sizeof(float *)*(col+1));
	for(i=0; i<(col+1); i++)
		Nx_imag[i] = (float *)malloc(sizeof(float)*(row+1));
	//-----------------------------------------------------------------------------------
	fourier_transform_Nx_Hy(row, col, dir, Nx_real, Nx_imag, hfieldy_real, hfieldy_imag);
	//-----------------------------------------------------------------------------------
	for(i=0; i<col; i++)
	{
		free(hfieldy_real[i]); free(hfieldy_imag[i]);
	}
		free(hfieldy_real); free(hfieldy_imag);
	//-------------------------------------------
	data_shifting_Nx(row, col, Nx_real, Nx_imag);
	//-------------------------------------------
	Nx_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Nx_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Nx_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Nx_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	//-------------------------------------------------------------------------
	Nx_data_shrink(row, col, Nx_real, Nx_imag, Nx_s_real, Nx_s_imag, Nfree, k);
	//-------------------------------------------------------------------------
	for(i=0; i<col+1; i++)
	{
		free(Nx_real[i]); free(Nx_imag[i]);
	}
		free(Nx_real); free(Nx_imag);

	Nt_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Nt_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Nt_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Nt_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1)); 
	Np_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Np_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Np_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Np_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Lt_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Lt_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Lt_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Lt_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Lp_s_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Lp_s_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	Lp_s_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		Lp_s_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	cosp = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		cosp[i] = (float *)malloc(sizeof(float)*(row_s+1));
	sinp = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		sinp[i] = (float *)malloc(sizeof(float)*(row_s+1));
	cost = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		cost[i] = (float *)malloc(sizeof(float)*(row_s+1));
    	ex_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		ex_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	ex_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		ex_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	ey_real = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		ey_real[i] = (float *)malloc(sizeof(float)*(row_s+1));
	ey_imag = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		ey_imag[i] = (float *)malloc(sizeof(float)*(row_s+1));
	printf("memory allocation ok! ...\n");

	///// old version of momentum space intensity distribution function 
	FT_data_writing(Nx_s_real, Nx_s_imag, Ny_s_real, Ny_s_imag, Lx_s_real, Lx_s_imag, Ly_s_real, Ly_s_imag,row_s, col_s,  Nfree, k);
	
	//--------------------------------------------------------------------------------------------------
	set_radiation_variable(Nfree, k, Eta, Nx_s_real, Nx_s_imag, Ny_s_real, Ny_s_imag, Lx_s_real, Lx_s_imag, Ly_s_real, Ly_s_imag, Nt_s_real, Nt_s_imag, Np_s_real, Np_s_imag, Lt_s_real, Lt_s_imag, Lp_s_real, Lp_s_imag, cosp, sinp, cost, ex_real, ex_imag, ey_real, ey_imag);
	//--------------------------------------------------------------------------------------------------
	for(i=0; i<col_s+1; i++)
	{
		free(Nx_s_real[i]); free(Nx_s_imag[i]); free(Ny_s_real[i]); free(Ny_s_imag[i]);
		free(Lx_s_real[i]); free(Lx_s_imag[i]); free(Ly_s_real[i]); free(Ly_s_imag[i]);
	}
		free(Nx_s_real); free(Nx_s_imag); free(Ny_s_real); free(Ny_s_imag);
		free(Lx_s_real); free(Lx_s_imag); free(Ly_s_real); free(Ly_s_imag);

	radint = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		radint[i] = (float *)malloc(sizeof(float)*(row_s+1));
	etheta = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		etheta[i] = (float *)malloc(sizeof(float)*(row_s+1));
	ephi = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		ephi[i] = (float *)malloc(sizeof(float)*(row_s+1));
	eintx = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		eintx[i] = (float *)malloc(sizeof(float)*(row_s+1));
	einty = (float **)malloc(sizeof(float *)*(col_s+1));
	for(i=0; i<(col_s+1); i++)
		einty[i] = (float *)malloc(sizeof(float)*(row_s+1));
	//--------------------------------------------------------------------------------------------------
	calc_radiation(Nfree, k, Eta, Nt_s_real, Nt_s_imag, Np_s_real, Np_s_imag, Lt_s_real, Lt_s_imag, Lp_s_real, Lp_s_imag, cosp, sinp, cost, ex_real, ex_imag, ey_real, ey_imag, radint, etheta, ephi, eintx, einty, &P_tot, &P_the, &P_phi);
	//--------------------------------------------------------------------------------------------------
	for(i=0; i<col_s+1; i++)
	{
		free(Nt_s_real[i]); free(Nt_s_imag[i]); free(Np_s_real[i]); free(Np_s_imag[i]);
		free(Lt_s_real[i]); free(Lt_s_imag[i]); free(Lp_s_real[i]); free(Lp_s_imag[i]);
	}
		free(Nt_s_real); free(Nt_s_imag); free(Np_s_real); free(Np_s_imag);
		free(Lt_s_real); free(Lt_s_imag); free(Lp_s_real); free(Lp_s_imag);

	//-----------------------------------------------------------------------------------------
	print_radiation(Nfree, k, Eta, radint, etheta, ephi, eintx, einty, &P_tot, &P_the, &P_phi);
	//-----------------------------------------------------------------------------------------

	//calc_polar(NA, Nfree);
	for(i=0; i<col_s+1; i++)
	{
		free(ex_real[i]); free(ey_real[i]); free(ex_imag[i]); free(ey_imag[i]);
	}
	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);
}

int FFT2D(struct COMPLEX **c,int nx,int ny,int dir)
{
   	int i,j;