output.c

program FDTD for the photonic crystal structure

			sprintf(name,"%07d",start);
			strcat(name,lastname);
			file=fopen(name,"w");

			for(k=k_range;k>=1;k--)
			{
			
				for(j=1;j<=j_range;j++)	fprintf(file,"%g ",(*field_avg)[j][k]/(float)(end-start));
				fprintf(file,"\n");
			}

			fclose(file);
		}
	}

	if(strcmp("y",plane)==0)
	{
		if(t == start)
		{
			///// *field_avg[k][i] /////
			(*field_avg) = (float **)malloc(sizeof(float *)*(k_range+1));
			for(k=0; k<k_range+1; k++)
				(*field_avg)[k] = (float *)malloc(sizeof(float)*(i_range+1));

			for(k=0; k<k_range+1; k++)
				for(i=0; i<i_range+1; i++)
					(*field_avg)[k][i] = 0.0;
		}

		if(t>start && t<=end)
		{
			j=floor(0.5+((value+ycenter)*lattice_y));
			for(k=k_range;k>=1;k--)
			{
				for(i=1;i<=i_range;i++)  
					(*field_avg)[k][i] = (*field_avg)[k][i] + grid_value(component,i,j,k); 
			}
		}

		if(t == end)
		{
			sprintf(name,"%07d",start);
			strcat(name,lastname);
			file=fopen(name,"w");

			for(k=k_range;k>=1;k--)
			{
			
				for(i=1;i<=i_range;i++)	fprintf(file,"%g ",(*field_avg)[k][i]/(float)(end-start));
				fprintf(file,"\n");
			}

			fclose(file);
		}
	}

	if(strcmp("z",plane)==0)
	{
		if(t == start)
		{
			///// *field_avg[i][j] /////
			(*field_avg) = (float **)malloc(sizeof(float *)*(i_range+1));
			for(i=0; i<i_range+1; i++)
				(*field_avg)[i] = (float *)malloc(sizeof(float)*(j_range+1));

			for(i=0; i<i_range+1; i++)
				for(j=0; j<j_range+1; j++)
					(*field_avg)[i][j] = 0.0;
		}

		if(t>start && t<=end)
		{
			k=non_uniform_z_to_i(value);
			for(j=j_range;j>=1;j--)
			{
				for(i=1;i<=i_range;i++)  
					(*field_avg)[i][j] = (*field_avg)[i][j] + grid_value(component,i,j,k); 
			}
		}

		if(t == end)
		{
			sprintf(name,"%07d",start);
			strcat(name,lastname);
			file=fopen(name,"w");

			for(j=j_range;j>=1;j--)
			{
			
				for(i=1;i<=i_range;i++)	fprintf(file,"%g ",(*field_avg)[i][j]/(float)(end-start));
				fprintf(file,"\n");
			}

			fclose(file);
		}
	}
}

void out_plane_time_average_projection(char *component,char *dirc, long start, long end, float ***field_avg, char *lastname, int k_shift)
{
	FILE *file;
	char name[20];
	int i,j,k;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	if(use_periodic_x == 1)
		i_range = isize;
	if(use_periodic_y == 1)
		j_range = jsize; 

	if(strcmp("+z",dirc)==0)
	{
		if(t == start)
		{
			///// *field_avg[i][j] /////
			(*field_avg) = (float **)malloc(sizeof(float *)*(i_range+1));
			for(i=0; i<i_range+1; i++)
				(*field_avg)[i] = (float *)malloc(sizeof(float)*(j_range+1));

			for(i=0; i<i_range+1; i++)
				for(j=0; j<j_range+1; j++)
					(*field_avg)[i][j] = 0.0;
		}

		if(t>start && t<=end)
		{
			for(j=j_range;j>=1;j--)
			{
				for(i=1;i<=i_range;i++)  
				{
					k= find_k_projection(i,j,k_range) + k_shift;
					(*field_avg)[i][j] = (*field_avg)[i][j] + grid_value(component,i,j,k); 
				}
			}
		}

		if(t == end)
		{
			sprintf(name,"%07d",start);
			strcat(name,lastname);
			file=fopen(name,"w");

			for(j=j_range;j>=1;j--)
			{
				for(i=1;i<=i_range;i++)
					fprintf(file,"%g ",(*field_avg)[i][j]/(float)(end-start));

				fprintf(file,"\n");
			}

			fclose(file);
		}
	}
}

void out_plane_projection(char *component,char *dirc,char *lastname, int k_shift)
{
	FILE *file;
	char name[20];
	int i,j,k;
	int i_range, j_range, k_range; 

	// for Hz_parity, normal cases
	i_range = isize-1;
	j_range = jsize-1;
	k_range = ksize-1; 
	
	if(use_periodic_x == 1)
		i_range = isize;
	if(use_periodic_y == 1)
		j_range = jsize; 
	
	sprintf(name,"%07d",t);
	strcat(name,lastname);
	file=fopen(name,"w");

	if(strcmp("+z",dirc)==0)
	{
		for(j=j_range;j>=1;j--)
		{
			for(i=1;i<=i_range;i++)
			{
				k= find_k_projection(i,j,k_range) + k_shift;
				fprintf(file,"%g ",grid_value(component,i,j,k));
			}	
			fprintf(file,"\n");
		}
	}
	fclose(file);
	printf("out %s...ok\n",component);
}

void out_several_points(char *component, float zposition, float xside, float yside, int pNx, int pNy, long ti, long tf, char *name)
{
	FILE *file;
	int i,j,k;
	int m,n; //for loop index

	//// xside, yside --> i_range, j_range
	
	if(ti<=t && t<tf)
	{
		k=non_uniform_z_to_i(zposition);  // fixed 

		file=fopen(name,"a");

		for(m=0; m<=pNx; m++)
		{
			for(n=0; n<=pNy; n++)
			{
				i=floor(0.5+((-0.5*xside+m*(xside/pNx)+xcenter)*lattice_x)); 
				j=floor(0.5+((-0.5*yside+n*(yside/pNx)+ycenter)*lattice_y));
				fprintf(file,"%g\t",grid_value(component,i,j,k));
			}
		}	
		fprintf(file,"\n");
		fclose(file);
	}

}

void out_point(char *component,float x,float y,float z,long ti,long tf,char *name)
{
	FILE *file;
	int i,j,k;

	if(ti<=t && t<tf)
	{
		i=floor(0.5+((x+xcenter)*lattice_x));
		j=floor(0.5+((y+ycenter)*lattice_y));
		k=non_uniform_z_to_i(z);

		if(t==ti) remove(name);
		file=fopen(name,"a");
		fprintf(file,"%g\n",grid_value(component,i,j,k));
		fclose(file);
	}
}

float grid_value(char *component,int i,int j,int k)
{
	if(strcmp(component,"Ex")==0) return efieldx(i,j,k);
	if(strcmp(component,"Ey")==0) return efieldy(i,j,k);
	if(strcmp(component,"Ez")==0) return efieldz(i,j,k);
	if(strcmp(component,"Hx")==0) return hfieldx(i,j,k);
	if(strcmp(component,"Hy")==0) return hfieldy(i,j,k);
	if(strcmp(component,"Hz")==0) return hfieldz(i,j,k);
	if(strcmp(component,"Jx")==0) return jfieldx(i,j,k);
	if(strcmp(component,"Jy")==0) return jfieldy(i,j,k);
	if(strcmp(component,"Jz")==0) return jfieldz(i,j,k);	
	if(strcmp(component,"J.E")==0) return efieldx(i,j,k)*jfieldx(i,j,k)+efieldy(i,j,k)*jfieldy(i,j,k)+efieldz(i,j,k)*jfieldz(i,j,k);
	if(strcmp(component,"E^2")==0) return pow(efieldx(i,j,k),2)+pow(efieldy(i,j,k),2)+pow(efieldz(i,j,k),2);
	if(strcmp(component,"H^2")==0) return pow(hfieldx(i,j,k),2)+pow(hfieldy(i,j,k),2)+pow(hfieldz(i,j,k),2);
	if(strcmp(component,"Sx")==0) return efieldy(i,j,k)*hfieldz(i,j,k)-efieldz(i,j,k)*hfieldy(i,j,k);
	if(strcmp(component,"Sy")==0) return efieldz(i,j,k)*hfieldx(i,j,k)-efieldx(i,j,k)*hfieldz(i,j,k);
	if(strcmp(component,"Sz")==0) return efieldx(i,j,k)*hfieldy(i,j,k)-efieldy(i,j,k)*hfieldx(i,j,k);
	if(strcmp(component,"LogE^2")==0) return log10(pow(efieldx(i,j,k),2)+pow(efieldy(i,j,k),2)+pow(efieldz(i,j,k),2));
	if(strcmp(component,"LogH^2")==0) return log10(pow(hfieldx(i,j,k),2)+pow(hfieldy(i,j,k),2)+pow(hfieldz(i,j,k),2));
	if(strcmp(component,"EM_Energy")==0) return eo*eps(i,j,k)*(pow(efieldx(i,j,k),2)+pow(efieldy(i,j,k),2)+pow(efieldz(i,j,k),2))+uo*ups*(pow(hfieldx(i,j,k),2)+pow(hfieldy(i,j,k),2)+pow(hfieldz(i,j,k),2));
	if(strcmp(component,"E_Energy")==0) return eo*eps(i,j,k)*(pow(efieldx(i,j,k),2)+pow(efieldy(i,j,k),2)+pow(efieldz(i,j,k),2));
	if(strcmp(component,"Ex2Ey2")==0) return pow(efieldx(i,j,k),2)+pow(efieldy(i,j,k),2);
	if(strcmp(component,"Ey2Ez2")==0) return pow(efieldy(i,j,k),2)+pow(efieldz(i,j,k),2);
	if(strcmp(component,"Ez2Ex2")==0) return pow(efieldz(i,j,k),2)+pow(efieldx(i,j,k),2);	
	if(strcmp(component,"LogEx2Ey2")==0) return log10(pow(efieldx(i,j,k),2)+pow(efieldy(i,j,k),2));
	if(strcmp(component,"LogEy2Ez2")==0) return log10(pow(efieldy(i,j,k),2)+pow(efieldz(i,j,k),2));
	if(strcmp(component,"LogEz2Ex2")==0) return log10(pow(efieldz(i,j,k),2)+pow(efieldx(i,j,k),2));
	if(strcmp(component,"Hx2Hy2")==0) return pow(hfieldx(i,j,k),2)+pow(hfieldy(i,j,k),2);
	if(strcmp(component,"Hy2Hz2")==0) return pow(hfieldy(i,j,k),2)+pow(hfieldz(i,j,k),2);
	if(strcmp(component,"Hz2Hx2")==0) return pow(hfieldz(i,j,k),2)+pow(hfieldx(i,j,k),2);	
	if(strcmp(component,"LogHx2Hy2")==0) return log10(pow(hfieldx(i,j,k),2)+pow(hfieldy(i,j,k),2));
	if(strcmp(component,"LogHy2Hz2")==0) return log10(pow(hfieldy(i,j,k),2)+pow(hfieldz(i,j,k),2));
	if(strcmp(component,"LogHz2Hx2")==0) return log10(pow(hfieldz(i,j,k),2)+pow(hfieldx(i,j,k),2));	
}

float grid_value_periodic_x(char *component,int i,int j,int k, int h)
{
	float E_x, E_y, E_z, J_x, J_y, J_z, H_x, H_y, H_z;

	E_x = efieldx(i,j,k)*cos(2*pi*h*wave_vector_x)-iefieldx(i,j,k)*sin(2*pi*h*wave_vector_x);
	E_y = efieldy(i,j,k)*cos(2*pi*h*wave_vector_x)-iefieldy(i,j,k)*sin(2*pi*h*wave_vector_x);
	E_z = efieldz(i,j,k)*cos(2*pi*h*wave_vector_x)-iefieldz(i,j,k)*sin(2*pi*h*wave_vector_x);
	J_x = jfieldx(i,j,k)*cos(2*pi*h*wave_vector_x)-ijfieldx(i,j,k)*sin(2*pi*h*wave_vector_x);
	J_y = jfieldy(i,j,k)*cos(2*pi*h*wave_vector_x)-ijfieldy(i,j,k)*sin(2*pi*h*wave_vector_x);
	J_z = jfieldz(i,j,k)*cos(2*pi*h*wave_vector_x)-ijfieldz(i,j,k)*sin(2*pi*h*wave_vector_x);
	H_x = hfieldx(i,j,k)*cos(2*pi*h*wave_vector_x)-ihfieldx(i,j,k)*sin(2*pi*h*wave_vector_x);
	H_y = hfieldy(i,j,k)*cos(2*pi*h*wave_vector_x)-ihfieldy(i,j,k)*sin(2*pi*h*wave_vector_x);
	H_z = hfieldz(i,j,k)*cos(2*pi*h*wave_vector_x)-ihfieldz(i,j,k)*sin(2*pi*h*wave_vector_x);

	if(strcmp(component,"Ex")==0) return (E_x);
	if(strcmp(component,"Ey")==0) return (E_y);
	if(strcmp(component,"Ez")==0) return (E_z);
	if(strcmp(component,"Hx")==0) return (H_x);
	if(strcmp(component,"Hy")==0) return (H_y);
	if(strcmp(component,"Hz")==0) return (H_z);
	if(strcmp(component,"Jx")==0) return (J_x);
	if(strcmp(component,"Jy")==0) return (J_y);
	if(strcmp(component,"Jz")==0) return (J_z);
	if(strcmp(component,"J.E")==0) return E_x*J_x+E_y*J_y+E_z*J_z;
	if(strcmp(component,"E^2")==0) return pow(E_x,2)+pow(E_y,2)+pow(E_z,2);
	if(strcmp(component,"H^2")==0) return pow(H_x,2)+pow(H_y,2)+pow(H_z,2);
	if(strcmp(component,"Sx")==0) return E_y*H_z-E_z*H_y;
	if(strcmp(component,"Sy")==0) return E_z*H_x-E_x*H_z;
	if(strcmp(component,"Sz")==0) return E_x*H_y-E_y*H_x;
	if(strcmp(component,"LogE^2")==0) return log10(pow(E_x,2)+pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"LogH^2")==0) return log10(pow(H_x,2)+pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"EM_Energy")==0) return eo*eps(i,j,k)*(pow(E_x,2)+pow(E_y,2)+pow(E_z,2))+uo*ups*(pow(H_x,2)+pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"E_Energy")==0) return eo*eps(i,j,k)*(pow(E_x,2)+pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"Ex2Ey2")==0) return pow(E_x,2)+pow(E_y,2);
	if(strcmp(component,"Ey2Ez2")==0) return pow(E_y,2)+pow(E_z,2);
	if(strcmp(component,"Ez2Ex2")==0) return pow(E_z,2)+pow(E_x,2);	
	if(strcmp(component,"LogEx2Ey2")==0) return log10(pow(E_x,2)+pow(E_y,2));
	if(strcmp(component,"LogEy2Ez2")==0) return log10(pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"LogEz2Ex2")==0) return log10(pow(E_z,2)+pow(E_x,2));
	if(strcmp(component,"Hx2Hy2")==0) return pow(H_x,2)+pow(H_y,2);
	if(strcmp(component,"Hy2Hz2")==0) return pow(H_y,2)+pow(H_z,2);
	if(strcmp(component,"Hz2Hx2")==0) return pow(H_z,2)+pow(H_x,2);	
	if(strcmp(component,"LogHx2Hy2")==0) return log10(pow(H_x,2)+pow(H_y,2));
	if(strcmp(component,"LogHy2Hz2")==0) return log10(pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"LogHz2Hx2")==0) return log10(pow(H_z,2)+pow(H_x,2));	
}

float grid_value_periodic_y(char *component,int i,int j,int k, int h)
{
	float E_x, E_y, E_z, J_x, J_y, J_z, H_x, H_y, H_z;

	E_x = efieldx(i,j,k)*cos(2*pi*h*wave_vector_y)-iefieldx(i,j,k)*sin(2*pi*h*wave_vector_y);
	E_y = efieldy(i,j,k)*cos(2*pi*h*wave_vector_y)-iefieldy(i,j,k)*sin(2*pi*h*wave_vector_y);
	E_z = efieldz(i,j,k)*cos(2*pi*h*wave_vector_y)-iefieldz(i,j,k)*sin(2*pi*h*wave_vector_y);
	J_x = jfieldx(i,j,k)*cos(2*pi*h*wave_vector_y)-ijfieldx(i,j,k)*sin(2*pi*h*wave_vector_y);
	J_y = jfieldy(i,j,k)*cos(2*pi*h*wave_vector_y)-ijfieldy(i,j,k)*sin(2*pi*h*wave_vector_y);
	J_z = jfieldz(i,j,k)*cos(2*pi*h*wave_vector_y)-ijfieldz(i,j,k)*sin(2*pi*h*wave_vector_y);
	H_x = hfieldx(i,j,k)*cos(2*pi*h*wave_vector_y)-ihfieldx(i,j,k)*sin(2*pi*h*wave_vector_y);
	H_y = hfieldy(i,j,k)*cos(2*pi*h*wave_vector_y)-ihfieldy(i,j,k)*sin(2*pi*h*wave_vector_y);
	H_z = hfieldz(i,j,k)*cos(2*pi*h*wave_vector_y)-ihfieldz(i,j,k)*sin(2*pi*h*wave_vector_y);

	if(strcmp(component,"Ex")==0) return (E_x);
	if(strcmp(component,"Ey")==0) return (E_y);
	if(strcmp(component,"Ez")==0) return (E_z);
	if(strcmp(component,"Hx")==0) return (H_x);
	if(strcmp(component,"Hy")==0) return (H_y);
	if(strcmp(component,"Hz")==0) return (H_z);
	if(strcmp(component,"Jx")==0) return (J_x);
	if(strcmp(component,"Jy")==0) return (J_y);
	if(strcmp(component,"Jz")==0) return (J_z);
	if(strcmp(component,"J.E")==0) return E_x*J_x+E_y*J_y+E_z*J_z;
	if(strcmp(component,"E^2")==0) return pow(E_x,2)+pow(E_y,2)+pow(E_z,2);
	if(strcmp(component,"H^2")==0) return pow(H_x,2)+pow(H_y,2)+pow(H_z,2);
	if(strcmp(component,"Sx")==0) return E_y*H_z-E_z*H_y;
	if(strcmp(component,"Sy")==0) return E_z*H_x-E_x*H_z;
	if(strcmp(component,"Sz")==0) return E_x*H_y-E_y*H_x;
	if(strcmp(component,"LogE^2")==0) return log10(pow(E_x,2)+pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"LogH^2")==0) return log10(pow(H_x,2)+pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"EM_Energy")==0) return eo*eps(i,j,k)*(pow(E_x,2)+pow(E_y,2)+pow(E_z,2))+uo*ups*(pow(H_x,2)+pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"E_Energy")==0) return eo*eps(i,j,k)*(pow(E_x,2)+pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"Ex2Ey2")==0) return pow(E_x,2)+pow(E_y,2);
	if(strcmp(component,"Ey2Ez2")==0) return pow(E_y,2)+pow(E_z,2);
	if(strcmp(component,"Ez2Ex2")==0) return pow(E_z,2)+pow(E_x,2);	
	if(strcmp(component,"LogEx2Ey2")==0) return log10(pow(E_x,2)+pow(E_y,2));
	if(strcmp(component,"LogEy2Ez2")==0) return log10(pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"LogEz2Ex2")==0) return log10(pow(E_z,2)+pow(E_x,2));
	if(strcmp(component,"Hx2Hy2")==0) return pow(H_x,2)+pow(H_y,2);
	if(strcmp(component,"Hy2Hz2")==0) return pow(H_y,2)+pow(H_z,2);
	if(strcmp(component,"Hz2Hx2")==0) return pow(H_z,2)+pow(H_x,2);	
	if(strcmp(component,"LogHx2Hy2")==0) return log10(pow(H_x,2)+pow(H_y,2));
	if(strcmp(component,"LogHy2Hz2")==0) return log10(pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"LogHz2Hx2")==0) return log10(pow(H_z,2)+pow(H_x,2));	
}

float grid_value_periodic_xy(char *component,int i,int j,int k, int h, int v)
{
	float E_x, E_y, E_z, J_x, J_y, J_z, H_x, H_y, H_z;

	E_x = efieldx(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-iefieldx(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	E_y = efieldy(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-iefieldy(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	E_z = efieldz(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-iefieldz(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	J_x = jfieldx(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-ijfieldx(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	J_y = jfieldy(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-ijfieldy(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	J_z = jfieldz(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-ijfieldz(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	H_x = hfieldx(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-ihfieldx(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	H_y = hfieldy(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-ihfieldy(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);
	H_z = hfieldz(i,j,k)*cos(2*pi*h*wave_vector_x+v*wave_vector_y)-ihfieldz(i,j,k)*sin(2*pi*h*wave_vector_x+v*wave_vector_y);

	if(strcmp(component,"Ex")==0) return (E_x);
	if(strcmp(component,"Ey")==0) return (E_y);
	if(strcmp(component,"Ez")==0) return (E_z);
	if(strcmp(component,"Hx")==0) return (H_x);
	if(strcmp(component,"Hy")==0) return (H_y);
	if(strcmp(component,"Hz")==0) return (H_z);
	if(strcmp(component,"Jx")==0) return (J_x);
	if(strcmp(component,"Jy")==0) return (J_y);
	if(strcmp(component,"Jz")==0) return (J_z);
	if(strcmp(component,"J.E")==0) return E_x*J_x+E_y*J_y+E_z*J_z;
	if(strcmp(component,"E^2")==0) return pow(E_x,2)+pow(E_y,2)+pow(E_z,2);
	if(strcmp(component,"H^2")==0) return pow(H_x,2)+pow(H_y,2)+pow(H_z,2);
	if(strcmp(component,"Sx")==0) return E_y*H_z-E_z*H_y;
	if(strcmp(component,"Sy")==0) return E_z*H_x-E_x*H_z;
	if(strcmp(component,"Sz")==0) return E_x*H_y-E_y*H_x;
	if(strcmp(component,"LogE^2")==0) return log10(pow(E_x,2)+pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"LogH^2")==0) return log10(pow(H_x,2)+pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"EM_Energy")==0) return eo*eps(i,j,k)*(pow(E_x,2)+pow(E_y,2)+pow(E_z,2))+uo*ups*(pow(H_x,2)+pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"E_Energy")==0) return eo*eps(i,j,k)*(pow(E_x,2)+pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"Ex2Ey2")==0) return pow(E_x,2)+pow(E_y,2);
	if(strcmp(component,"Ey2Ez2")==0) return pow(E_y,2)+pow(E_z,2);
	if(strcmp(component,"Ez2Ex2")==0) return pow(E_z,2)+pow(E_x,2);	
	if(strcmp(component,"LogEx2Ey2")==0) return log10(pow(E_x,2)+pow(E_y,2));
	if(strcmp(component,"LogEy2Ez2")==0) return log10(pow(E_y,2)+pow(E_z,2));
	if(strcmp(component,"LogEz2Ex2")==0) return log10(pow(E_z,2)+pow(E_x,2));
	if(strcmp(component,"Hx2Hy2")==0) return pow(H_x,2)+pow(H_y,2);
	if(strcmp(component,"Hy2Hz2")==0) return pow(H_y,2)+pow(H_z,2);
	if(strcmp(component,"Hz2Hx2")==0) return pow(H_z,2)+pow(H_x,2);	
	if(strcmp(component,"LogHx2Hy2")==0) return log10(pow(H_x,2)+pow(H_y,2));
	if(strcmp(component,"LogHy2Hz2")==0) return log10(pow(H_y,2)+pow(H_z,2));
	if(strcmp(component,"LogHz2Hx2")==0) return log10(pow(H_z,2)+pow(H_x,2));	
}

float efieldx(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(Ex[i][j][k]+Ex[i-1][j][k])/2;
}

float efieldy(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(Ey[i][j][k]+Ey[i][j-1][k])/2;
}

float efieldz(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*zparity;}

	return p*(Ez[i][j][k]+Ez[i][j][k-1])/2;
}

float jfieldx(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*xparity;}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*(-yparity);}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}

	return p*(Jx[i][j][k]+Jx[i-1][j][k])/2;
}

float jfieldy(int i,int j,int k)
{
	int p=1;

	if(pisize<=i && (xparity==1 || xparity==-1)) {i=isize-i;p=p*(-xparity);}
	if(pjsize<=j && (yparity==1 || yparity==-1)) {j=jsize-j;p=p*yparity;}
	if(pksize<=k && (zparity==1 || zparity==-1)) {k=ksize-k;p=p*(-zparity);}