/*fd2d_3.4.c. 2D TM simulation of a plane wave source impinging on a dielectric cylinder
              Analysis using Fourier Transforms*/

# include <math.h>
# include <stdlib.h>
# include <stdio.h>

#define IE 60
#define JE 60
#define NFREQS 3

void main()
{
	double amp[IE][JE]={0};
	double amp_in[NFREQS],phase_in[NFREQS];
	double radius=10,epsilon=30,sigma=0.3;
	double ga[IE][JE],gb[IE][JE],ez[IE][JE],hx[IE][JE],hy[IE][JE];
	double freq[NFREQS],arg[NFREQS];
	int n,i,j,nsteps,ia,ib,ic,ja,jb,jc,m;
	double ddx,dt,T,pi,epsz;
	double xdist,ydist,dist;
	double curl_e;
	double t0,spread,pulse;
	double gi2[IE],gi3[IE];
	double gj2[JE],gj3[JE];
	double fi1[IE],fi2[IE],fi3[IE];
	double fj1[IE],fj2[JE],fj3[JE];
//	double ihx[IE][JE],ihy[IE][JE];
	FILE *fp;
	ic=IE/2;
	jc=JE/2;
	ia=7;
	ib=IE-ia-1;
	ja=7;
	jb=JE-ja-1;
	ddx=0.01;                         /*Cell size*/
	dt=ddx/6e8;                       /*Time steps*/            
	epsz=8.8e-12;
	pi=3.14159;
    double real_in[NFREQS]={0},imag_in[NFREQS]={0}; 
	double ihx[IE][JE]={0},ihy[IE][JE]={0},iz[IE][JE]={0},dz[IE][JE]={0};
	double real_pt[NFREQS][IE][JE]={0},imag_pt[NFREQS][IE][JE]={0};
	double ezinc[IE],hxinc[JE];
	for(j=0;j<JE;j++)
	{
		ezinc[j]=0.0;
	}
    for(j=0;j<JE;j++)
	{
        hxinc[j]=0.0;
	}
//	printf("%2d,%2d",ic,jc);

	/*Initialize the arrays*/
	for(j=0;j<JE;j++) 
	{
//	printf("%2d",j);
	for(i=0;i<IE;i++) 
	{
		dz[i][j]=0;
		ez[i][j]=0;
		hx[i][j]=0;
		hy[i][j]=0;
		ga[i][j]=1.0;
		gb[i][j]=0.0;
	//	printf("%5.2f",ga[i][j]);
	}
//	printf("\n");
        }


	/*Calculate the PML parameters*/
	for(i=0;i<IE;i++) {
		gi2[i]=1.0;
		gi3[i]=1.0;
		fi1[i]=0.0;
		fi2[i]=1.0;
		fi3[i]=1.0;
	}

	for(j=0;j<JE;j++) {
		gj2[j]=1.0;
		gj3[j]=1.0;
		fj1[j]=0.0;
		fj2[j]=1.0;
		fj3[j]=1.0;
	}

	/* Parameters for the Fourier Transforms */
	freq[0]=50.e6;
	freq[1]=300.e6;
	freq[2]=700.e6;

	for(n=0;n<NFREQS;n++)
	{arg[n]=2*pi*freq[n]*dt;
	 printf("%d %6.2f %7.5 \n", n,freq[n]*1e-6,arg[n]);
	}

	/*Specify the dielectric cylinder*/
	printf("Cylinder radius (cells),epsilon,sigma-->");
//	scanf("%f",&radius);
//	scanf("%f",&epsilon);
//	scanf("%f",&sigma);
	printf("Radius= %5.2f Eps= %6.2f Sigma= %6.2f \n",radius,epsilon,sigma);

	for(j=ja;j<jb;j++) {
		for(i=ia;i<ib;i++) {
			xdist = (ic-i);
			ydist = (jc-j);
			dist = sqrt(pow(xdist,2.0)+pow(ydist,2.0));
			if(dist<=radius) {
				ga[i][j]=1.0/(epsilon+(sigma*dt/epsz));
				gb[i][j]=sigma*dt/epsz;
			}
		}
	}

	printf("\n");

/*	printf("Ga \n");
	for(j=ja;j<jb;j++)
	{   
		printf("%d \n",j);
		for(i=ia;i<ib;i++)
		{
			printf("%5.2f ",ga[i][j]);
		}
		
		printf("\n");
	}
	
	printf("\n");

	printf("Gb \n");
	for(j=ja;j<jb;j++)
	{
		printf("%d \n",j);
    	for(i=ia;i<ib;i++) 
		{
			printf("%5.2f ",gb[i][j]);
		}
		
		printf("\n");
	}

  */

    t0=25.0;
	spread=8.0;
	T=0;
	nsteps=1;

  // 	while(nsteps>0) {
		printf("nsteps-->");
		scanf("%d",&nsteps);
//		printf("%d \n",nsteps);

		for(n=1;n<=nsteps;n++){
			T=T+1;

	/* Start of the Main FDTD loop */
			
			/*incident loop 1D */
			for(j=1;j<JE;j++) {
				ezinc[j]=ezinc[j]+0.5*(hxinc[j-1]-hxinc[j]);
			}

   /*Fourier Transform of the incident field*/
	for(m=0;m<NFREQS;m++)
	{ real_in[m]=real_in[m]+cos(arg[m]*T)*ezinc[ja-1];
	  imag_in[m]=imag_in[m]+sin(arg[m]*T)*ezinc[ja-1];
	}


    /* Calculate the Dz filed*/
	for(j=1;j<JE;j++) {
		for(i=1;i<IE;i++){
			dz[i][j]=gi3[i]*gj3[j]*dz[i][j]+gi2[i]*gj2[j]*0.5*(hy[i][j]-hy[i-1][j]-hx[i][j]+hx[i][j-1]);		    				
		}
		}
    	
	printf("dz2020= %5.2f hy2020= %5.2f hy1920= %5.2f hy2019= %5.2f",dz[20][20],hy[20][20],hy[19][20],hy[20][19]);
				printf("\n");	


	/*sinusoidal Source*/
    pulse=exp(-0.5*(pow((t0-T)/spread,2.0)));
	ezinc[3]=pulse;
    printf("%3.0f %6.2f \n",T,ezinc[3]);

	/*Incident Dz values*/
		for(i=ia;i<=ib;i++) {
			dz[i][ja]=dz[i][ja]+0.5*hxinc[ja-1];
			dz[i][jb]=dz[i][jb]-0.5*hxinc[jb];
		}

     /*Calculate the Ez field*/
        for(j=1;j<JE-1;j++) {
			for(i=1;i<IE-1;i++) {
				ez[i][j]=ga[i][j]*(dz[i][j]-iz[i][j]);
				iz[i][j]=iz[i][j]+gb[i][j]*ez[i][j];
			}
		}
	
		
		
	/*Calculate the Fourier transform of Ez*/
	for(j=0;j<JE;j++) {   
		for(i=0;i<JE;i++){   
			for(m=0;m<NFREQS;m++){   
				real_pt[m][i][j]=real_pt[m][i][j]+cos(arg[m]*T)*ez[i][j];
				imag_pt[m][i][j]=imag_pt[m][i][j]+sin(arg[m]*T)*ez[i][j];
			}
		}
	}


	
		/*incident loop 1D*/
		for(j=0;j<JE-1;j++) {
			hxinc[j]=hxinc[j]+0.5*(ezinc[j]-ezinc[j+1]);
		}       



		/*Calculate the Hx field*/
		for(j=0;j<JE-1;j++) {
			for(i=0;i<IE-1;i++) {
				curl_e=ez[i][j]-ez[i][j+1];
				ihx[i][j]=ihx[i][j]+curl_e;
				hx[i][j]=fj3[j]*hx[i][j]+fj2[j]*0.5*(curl_e+fi1[i]*ihx[i][j]);
			}
		}

		/*Incident hx values*/
		for(i=ia;i<=ib;i++) {
			hx[i][ja-1]=hx[i][ja-1]+0.5*ezinc[ja];
			hx[i][jb]=hx[i][jb]-0.5*ezinc[jb];
		}

		/*Calculate the Hy field*/
		for(j=0;j<JE-1;j++) {
			for(i=0;i<IE-1;i++) {
				curl_e=ez[i+1][j]-ez[i][j];
				ihy[i][j]=ihy[i][j]+curl_e;
				hy[i][j]=fi3[i]*hy[i][j]+fi2[i]*0.5*(curl_e+fj1[j]*ihy[i][j]);
			}
		}

		/*Incident hy values*/
		for(j=ja;j<=jb;j++) {
			hy[ia-1][j]=hy[ia-1][j]-0.5*ezinc[j];
			hy[ib][j]=hy[ib][j]+0.5*ezinc[j];
		}
		
		
		
		}

			/*End of the main FDTD loop*/

			
     /*Write the E field out to a file "Ez"*/
   fp=fopen("Ez.dat","w");
	for(j=0;j<JE;j++) {
		for(i=0;i<IE;i++) {
			fprintf(fp,"%6.3f ",ez[i][j]);
		}
		fprintf(fp,"\n");
	}
	fclose(fp);   
    printf("T= %6.0f \n",T);   

	/*Calculate the Fouier amplitude and phase of the incident pulse*/
/*	for(m=0;m<NFREQS;m++) {
		amp_in[m]=sqrt(pow(real_in[m],2.0)+pow(imag_in[m],2.0));
		phase_in[m]=atan2(imag_in[m],real_in[m]);
		printf("%d input Pulse: %8.4f %8.4f %8.4f %7.2f\n",m,real_in[m],imag_in[m],amp_in[m],(180.0/pi)*phase_in[m]);
	}
*/
    /*Calculate the Fouier amplitude and phase of the total field*/
/*	for (m=0;m<NFREQS;m++)
	{
		if(m==0) fp=fopen("amp1","w");
		else if(m==1) fp=fopen("amp2","w");
		else if (m==2) fp=fopen("amp3","w");
		{
			printf("%2d %7.2f MHz\n",m,freq[m]*1e-6);
			for(j=ja;j<=jb;j++) 
			{
				if(ga[ic][j]<1.00)
				{    amp[ic][j]=(1.0/amp_in[m])*sqrt(pow(real_pt[m][ic][j],2.0)+pow(imag_pt[m][ic][j],2.0));
				     printf("%2d %9.4f \n",jc-j,amp[ic][j]);
				     fprintf(fp,"%d %9.4f \n", j,amp[ic][j]);
				}
		
			}
		}
		fclose(fp);
	}*/
}