fd2d_3.3.cpp

2D TM program with plane wave source fd2d_3.3.c

//fd2d_3.3.c   2D TM program with plane wave source

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

#define IE 60
#define JE 60

main()
{
	float ez_inc[JE], hx_inc[JE];
	float ez_inc_low_m1, ez_inc_low_m2;
	float ez_inc_high_m1, ez_inc_high_m2;
	int ia, ib, ja, jb;

	float ga[IE][JE], dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE];
	int l, n, i, j, ic, jc, nsteps, npml;
	float ddx, dt, T, epsz, pi, epsilon, sigma, eaf;
	float xn, xxn, xnum, xd, curl_e;
	float t0, spread, pulse;
	float gi2[IE], gi3[IE];
	float gj2[JE], gj3[JE];
	float fi1[IE],fi2[IE],fi3[IE];
	float fj1[JE],fj2[JE],fj3[JE];
	float ihx[IE][JE], ihy[IE][JE];

	FILE *fp;

	ic=IE/2-15;
	jc=JE/2-15;

	ia=7;               //Total/scattered field boundaries
	ib=IE-ia-1;
	ja=7;
	jb=JE-ja-1;

	ddx=0.01;     //cell size
	dt=ddx/6e8;   //time step
	epsz=8.8e-12;
	pi=3.14159;

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

	for (j=0; j<JE; j++){
		ez_inc[j]=0, hx_inc[j]=0;
	}
	ez_inc_low_m1=0, ez_inc_low_m2=0;
	ez_inc_high_m1=0, ez_inc_high_m2=0;



	/*----------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;
	}

	printf("Number of PML cells --> ");
	scanf("%d", &npml);

	for(i=0;i<=npml;i++){
		xnum=npml-i;
		xd=npml;
		xxn=xnum/xd;
		xn=0.33*pow(xxn,3.0);
		printf(" %d   %7.4f   %7.4f  \n", i, xxn, xn);
		gi2[i]=1.0/(1.0+xn);
		gi2[IE-1-i]=1.0/(1.0+xn);
		gi3[i]=(1.0-xn)/(1.0+xn);
		gi3[IE-1-i]=(1.0-xn)/(1.0+xn);
		xxn=(xnum-0.5)/xd;   //orig
		//xxn=(xnum)/xd;

		xn=0.25*pow(xxn,3.0);
		fi1[i]=xn;
		fi1[IE-2-i]=xn;
		fi2[i]=1.0/(1.0+xn);
		fi2[IE-2-i]=1.0/(1.0+xn);
		fi3[i]=(1.0-xn)/(1.0+xn);
		fi3[IE-2-i]=(1.0-xn)/(1.0+xn);
	}

	for(j=0;j<=npml;j++){
		xnum=npml-j;
		xd=npml;
		xxn=xnum/xd;
		xn=0.33*pow(xxn,3.0);
		printf("%d  %7.4f  %7.4f  \n", j, xxn, xn);
		gj2[j]=1.0/(1.0+xn);
		gj2[JE-1-j]=1.0/(1.0+xn);
		gj3[j]=(1.0-xn)/(1.0+xn);
		gj3[JE-1-j]=(1.0-xn)/(1.0+xn);
		xxn=(xnum-0.5)/xd;  //orig
		//xxn=(xnum)/xd;

		xn=0.25*pow(xxn,3.0);
		fj1[j]=xn;
		fj1[JE-2-j]=xn;
		fj2[j]=1.0/(1.0+xn);
		fj2[JE-2-j]=1.0/(1.0+xn);
		fj3[j]=(1.0-xn)/(1.0+xn);
		fj3[JE-2-j]=(1.0-xn)/(1.0+xn);
	}




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

	while (nsteps>0){
		printf("nsteps --> "); 
		scanf("%d", &nsteps);
		for (n=1; n<=nsteps; n++){
			T=T+1;

			/* start of the Main loop */
			for (j=1; j<JE-1;j++){
				ez_inc[j]=ez_inc[j]+0.5*(hx_inc[j-1]-hx_inc[j]);
			}

			//ABC for the incident buffer
			
			ez_inc[0]    =ez_inc_low_m2;
			ez_inc_low_m2=ez_inc_low_m1;
			ez_inc_low_m1=ez_inc[1];
			//ez_inc_low_m2=ez_inc_low_m1;
			//ez_inc[0]    =ez_inc_low_m2;

			ez_inc[JE-1]    =ez_inc_high_m2;
			ez_inc_high_m2=ez_inc_high_m1;
			ez_inc_high_m1=ez_inc[JE-2];
			//ez_inc_high_m2=ez_inc_high_m1;
			//ez_inc[JE-1]    =ez_inc_high_m2;
			

			//calculate the Dz field
			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]);
				}
			}

			//source
			pulse=exp(-0.5*( pow((t0-T)/spread,2.0) ));
			ez_inc[3]=pulse;

			//Incident Dz value
			for (i=ia;i<=ib;i++){
				dz[i][ja]=dz[i][ja]+0.5*hx_inc[ja-1];
				dz[i][jb]=dz[i][jb]-0.5*hx_inc[jb];
			}
			
//(1)
			//calculate the Ez field
			for(j=0;j<JE;j++){
				//for(i=0;i<IE;i++){
					for(i=0;i<IE;i++){
					ez[i][j]=ga[i][j]*dz[i][j];
				}
			}
//(2)
			/*   set the Ez edges to 0,  as part of the PML  */
			for(j=0;j<JE-1;j++){
				ez[0][j]=0.0;
				ez[IE-1][j]=0.0;
			}
			for(i=0;i<IE-1;i++){
				ez[i][0]=0.0;
				ez[i][JE-1]=0.0;
			}

			for(j=0;j<JE;j++){
				hx_inc[j]=hx_inc[j]+0.5*(ez_inc[j]-ez_inc[j+1]);
			}

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

			//incident Hx value
			for (i=ia;i<=ib;i++){
				hx[i][ja-1] =hx[i][ja-1] + 0.5*ez_inc[ja];
				hx[i][jb]   =hx[i][jb]   - 0.5*ez_inc[jb];
			}

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

			//incident Hy value
			for (j=ja;j<=jb;j++){
				hy[ia-1][j] =hy[ia-1][j] - 0.5*ez_inc[j];
				hy[ib][j]   =hy[ib][j]   + 0.5*ez_inc[j];
			}

		}
			

	/*-----------End of the main FDTD loop-------------*/

	for(j=1;j<jc;j++){
		printf("%2d  ", j);
		for(i=1;i<ic;i++){
			printf("%5.2f ",ez[2*i][2*j]);
		}
		printf("\n");
	}

	//printf("T=%5.0f \n",T);

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

}