fd2d_3_l_01.cpp

采用fdtd计算二维电磁场

/* fd2d-3.l.c. 2D TM program */
# include <math.h>
# include <stdlib.h>
# include <stdio.h>

#define IE 60
#define JE 60

main()
{
	double ga[IE][JE] , dz[IE][JE] , ez[IE][JE] , hx[IE][JE] , hy [IE][JE];
	int n, i, j, ic, jc, nsteps;	//l, 
	double ddx, dt, T, epsz, pi;	// epsilon, eaf;, sigma
	double t0, spread, pulse;
	FILE *fp;

	ic = IE/2;
	jc = JE/2;
	ddx = 0.01;				/* Cell size */
	dt = ddx / 2 * 3.0e8;	/* Time steps */
	epsz = 8.85419e-12;
	pi = 3.1415926535;
	fp = fopen( "Ez", "w" );
	
	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;
			ga[i][j] = 1.0;
			printf( "%5.2f ", ga[i][j]);
		}
		printf(" \n");
	}

	t0 = 20.0;
	spread = 6.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 - - - - * /
			/* Calculate the Dz f i e l d */
			for( j = 1 ; j < IE; j++ )
			{
				for( i = 1 ; i < IE; i++ )
				{
					dz[i][j] = dz[i][j] + 0.5 * ( hy[i][j] - hy[i-1][j] - hx[i][j] + hx[i][j-1]);
				}
			}
			/* Put a Gaussian pulse in the middle */
			pulse = exp( -0.5 * ( pow ( ( t0 - T ) / spread , 2.0 ) ) );
			dz[ic][jc] = pulse;
			/* Calculate the Ez field */
			for( j = 1 ; j < JE; j++ )
			{
				for( i = 1 ; i < IE ; i++ )
				{
					ez[i][j] = ga[i][j] * dz[i][j];
				}
			}
			/* Calculate the Hx field */
			for( j=0; j < JE-1; j++ )
			{
				for( i = 0 ; i < IE-1; i++)
				{
					hx[i][j] = hx[i][j] + 0.5 * ( ez[i][j] - ez[i][j+1]);
				}
			}
			/* Calculate the Hy field */
			for ( j=0; j < JE-1; j++ )
			{
				for( i=0; i <= IE-1; i++ )
				{
					hy[i][j] = hy[i][j] + 0.5 * ( ez[i+1][j] - ez[i][j]);
				}
			}
			//Print the computational result of E field
			/* Write the E field out to a file "Ez" */
			fprintf( fp, "Time step = %5.0f \n[", T );
			for( j=0; j < JE; j++ )
			{
				for( i=0; i < IE; i++ )
				{
					fprintf( fp, "%6.3f   ", ez[i][j] );
				}
				fprintf( fp, ";\n") ;
			}
			fprintf( fp, "]\n");
			fflush( fp );
		}
		/* - - - - 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 );
//		fprintf( fp, "[" ) ;
//		fprintf( fp, "]" ) ;

		fclose(fp);
	}
	return 0;
}