calculation.c

利用不均一分割的FDTD解析程序

/*
 *      Calculation of Electric Field Componets
 *           and Magnetic Field Components
 */

void calculation(double vz)
{

	int i,j,k;
	double aax,aay,aaz,bbx,bby,bbz,r1,r2,r3,r4;
	double delta_x,delta_y,delta_z;
	double ux31,ux41,uy31,uy41,uz31,uz41;
	double hx,hy,hz;

    /******   Incident Pulse (Gap Excitation)   ******/
	ezi[Fx][Fy][Fz] = -vz/dz_sub;

	/******  Compting Magnetic Field Components  ******/
	for(i=0;i<=xs;i++){
		for(j=0;j<=(ys-1);j++){
			delta_y = dy;

			if((j >= (Fy-y_sub_length))
			    &&(j <= (Fy+y_sub_length-1))){
				delta_y = dy_sub;
			}

			for(k=0;k<=(zs-1);k++){
				delta_z = dz;

				if(k <= (Fz+z_sub_length-1)){
					delta_z = dz_sub;
				}

				aax = (eyi[i][j][k+1]-eyi[i][j][k])/delta_z
				    -(ezi[i][j+1][k]-ezi[i][j][k])/delta_y;
				hxi[i][j][k] = hxi[i][j][k]+rh1*aax;

				if(k == (Fz+z_sub_length)){
					hx_z_absorb[i][j][1]=hxi[i][j][k-1]*Inter11
					    +hxi[i][j][k]*Inter12;
				}

				if(j == (Fy-y_sub_length)){
					hx_y_absorb[i][0][k]= hxi[i][j-2][k]*Inter21
					    +hxi[i][j-1][k]*Inter22;
				}
				if(j == (Fy+y_sub_length)){
					hx_y_absorb[i][1][k]=hxi[i][j-1][k]*Inter11
					    +hxi[i][j][k]*Inter12;
				}
			}
		}
	}


	for(i=0;i<=(xs-1);i++){
		delta_x = dx;

		if((i >= (Fx-x_sub_length))
		    &&(i <= (Fx+x_sub_length-1))){
			delta_x = dx_sub;
		}

		for(j=0;j<=ys;j++){
			for(k=0;k<=(zs-1);k++){
				delta_z = dz;

				if(k <= (Fz+z_sub_length-1)){
					delta_z = dz_sub;
				}

				aay = (ezi[i+1][j][k]-ezi[i][j][k])/delta_x
				    -(exi[i][j][k+1]-exi[i][j][k])/delta_z;
				hyi[i][j][k] = hyi[i][j][k]+rh1*aay;

				if(k == (Fz+z_sub_length)){
					hy_z_absorb[i][j][1] = hyi[i][j][k-1]*Inter11
					    +hyi[i][j][k]*Inter12;
				}
				if(i == (Fx-x_sub_length)){
					hy_x_absorb[0][j][k] = hyi[i-2][j][k]*Inter21
					    +hyi[i-1][j][k]*Inter22;
				}
				if(i == (Fx+x_sub_length)){
					hy_x_absorb[1][j][k] = hyi[i-1][j][k]*Inter11
					    +hyi[i][j][k]*Inter12;
				}
			}
		}
	}

	for(i=0;i<=(xs-1);i++){
		delta_x = dx;

		if((i >= (Fx-x_sub_length))
		    &&(i <= (Fx+x_sub_length-1))){
			delta_x = dx_sub;
		}

		for(j=0;j<=(ys-1);j++){
			delta_y = dy;

			if((j >= (Fy-y_sub_length))
			    &&(j <= (Fy+y_sub_length-1))){
				delta_y = dy_sub;
			}

			for(k=0;k<=zs;k++){
				aaz = (exi[i][j+1][k]-exi[i][j][k])/delta_y
				    -(eyi[i+1][j][k]-eyi[i][j][k])/delta_x;
				hzi[i][j][k] = hzi[i][j][k]+rh1*aaz;

				if(i == (Fx-x_sub_length)){
					hz_x_absorb[0][j][k] = hzi[i-2][j][k]*Inter21
					    +hzi[i-1][j][k]*Inter22;
				}
				if(i == (Fx+x_sub_length)){
					hz_x_absorb[1][j][k] = hzi[i-1][j][k]*Inter11
					    +hzi[i][j][k]*Inter12;
				}
				if(j == (Fy-y_sub_length)){
					hz_y_absorb[i][0][k] = hzi[i][j-2][k]*Inter21
					    +hzi[i][j-1][k]*Inter22;
				}
				if(j == (Fy+y_sub_length)){
					hz_y_absorb[i][1][k] = hzi[i][j-1][k]*Inter11
					    +hzi[i][j][k]*Inter12;
				}
			}
		}
	}


   /*
	*    Saving Electric Field Components Temporary
	*    for Absorbing Boundary Condition   
	*/

	for(i=0;i<=(xs-1);i++){
		for(k=0;k<=zs;k++){
			exiya2[i][0][k] = exiya1[i][0][k];
			exiya2[i][1][k] = exiya1[i][1][k];
			exiyb2[i][ys][k] = exiyb1[i][ys][k];
			exiyb2[i][ys-1][k] = exiyb1[i][ys-1][k];
			exiya1[i][0][k] = exi[i][0][k];
			exiya1[i][1][k] = exi[i][1][k];
			exiyb1[i][ys][k] = exi[i][ys][k];
			exiyb1[i][ys-1][k] = exi[i][ys-1][k];
		}
	}

	for(i=0;i<=(xs-1);i++){
		for(j=0;j<=ys;j++){
			exiza2[i][j][0] = exiza1[i][j][0];
			exiza2[i][j][1] = exiza1[i][j][1];
			exizb2[i][j][zs] = exizb1[i][j][zs];
			exizb2[i][j][zs-1] = exizb1[i][j][zs-1];
			exiza1[i][j][0] = exi[i][j][0];
			exiza1[i][j][1] = exi[i][j][1];
			exizb1[i][j][zs] = exi[i][j][zs];
			exizb1[i][j][zs-1] = exi[i][j][zs-1];
		}
	}

	for(i=0;i<=xs;i++){
		for(j=0;j<=(ys-1);j++){
			eyiza2[i][j][0] = eyiza1[i][j][0];
			eyiza2[i][j][1] = eyiza1[i][j][1];
			eyizb2[i][j][zs] = eyizb1[i][j][zs];
			eyizb2[i][j][zs-1] = eyizb1[i][j][zs-1];
			eyiza1[i][j][0] = eyi[i][j][0];
			eyiza1[i][j][1] = eyi[i][j][1];
			eyizb1[i][j][zs] = eyi[i][j][zs];
			eyizb1[i][j][zs-1] = eyi[i][j][zs-1];
		}
	}

	for(j=0;j<=(ys-1);j++){
		for(k=0;k<=zs;k++){
			eyixa2[0][j][k] = eyixa1[0][j][k];
			eyixa2[1][j][k] = eyixa1[1][j][k];
			eyixb2[xs][j][k] = eyixb1[xs][j][k];
			eyixb2[xs-1][j][k] = eyixb1[xs-1][j][k];
			eyixa1[0][j][k] = eyi[0][j][k];
			eyixa1[1][j][k] = eyi[1][j][k];
			eyixb1[xs][j][k] = eyi[xs][j][k];
			eyixb1[xs-1][j][k] = eyi[xs-1][j][k];
		}
	}

	for(j=0;j<=ys;j++){
		for(k=0;k<=(zs-1);k++){
			ezixa2[0][j][k] = ezixa1[0][j][k];
			ezixa2[1][j][k] = ezixa1[1][j][k];
			ezixb2[xs][j][k] = ezixb1[xs][j][k];
			ezixb2[xs-1][j][k] = ezixb1[xs-1][j][k];
			ezixa1[0][j][k] = ezi[0][j][k];
			ezixa1[1][j][k] = ezi[1][j][k];
			ezixb1[xs][j][k] = ezi[xs][j][k];
			ezixb1[xs-1][j][k] = ezi[xs-1][j][k];
		}
	}

	for(i=0;i<=xs;i++){
		for(k=0;k<=(zs-1);k++){
			eziya2[i][0][k] = eziya1[i][0][k];
			eziya2[i][1][k] = eziya1[i][1][k];
			eziyb2[i][ys][k] = eziyb1[i][ys][k];
			eziyb2[i][ys-1][k] = eziyb1[i][ys-1][k];
			eziya1[i][0][k] = ezi[i][0][k];
			eziya1[i][1][k] = ezi[i][1][k];
			eziyb1[i][ys][k] = ezi[i][ys][k];
			eziyb1[i][ys-1][k] = ezi[i][ys-1][k];
		}
	}


	/*****   Computing Electric Field Components   *****/

	for(i=0;i<=(xs-1);i++){
		for(j=1;j<=(ys-1);j++){
			delta_y = dy;

			if((j >= (Fy-y_sub_length))
			    &&(j <= (Fy+y_sub_length))){
				delta_y = dy_sub;
			}

			for(k=1;k<=(zs-1);k++){
				delta_z = dz;
				hy = hyi[i][j][k]-hyi[i][j][k-1];

				if(k <= (Fz+z_sub_length)){
					delta_z = dz_sub;

					if(k == (Fz+z_sub_length)){
						hy=hy_z_absorb[i][j][1]-hyi[i][j][k-1];
					}
				}

				hz = hzi[i][j][k]-hzi[i][j-1][k];
				if(j == (Fy-y_sub_length)){
					hz = hzi[i][j][k]-hz_y_absorb[i][0][k];
				}
				if(j == (Fy+y_sub_length)){
					hz = hz_y_absorb[i][1][k]-hzi[i][j-1][k];
				}

				bbx = hz/delta_y-hy/delta_z;
				exi[i][j][k] = exi[i][j][k]+rh2*bbx;
			}
		}
	}

	for(i=1;i<=(xs-1);i++){
		delta_x = dx;

		if((i >= (Fx-x_sub_length))
		    &&(i <= (Fx+x_sub_length))){
			delta_x = dx_sub;
		}

		for(j=0;j<=(ys-1);j++){
			for(k=1;k<=(zs-1);k++){
				delta_z = dz;

				hx = hxi[i][j][k]-hxi[i][j][k-1];
				if(k <= (Fz+z_sub_length)){
					delta_z = dz_sub;
					if(k == (Fz+z_sub_length)){
						hx = hx_z_absorb[i][j][1]-hxi[i][j][k-1];
					}
				}

				hz = hzi[i][j][k]-hzi[i-1][j][k];
				if(i == (Fx-x_sub_length)){
					hz = hzi[i][j][k]-hz_x_absorb[0][j][k];
				}
				if(i == (Fx+x_sub_length)){
					hz = hz_x_absorb[1][j][k]-hzi[i-1][j][k];
				}

				bby = hx/delta_z-hz/delta_x;
				eyi[i][j][k] = eyi[i][j][k]+rh2*bby;
			}
		}
	}


	for(i=1;i<=(xs-1);i++){
		delta_x = dx;

		if((i >= (Fx-x_sub_length))
		    &&(i <= (Fx+x_sub_length))){
			delta_x = dx_sub;
		}

		for(j=1;j<=(ys-1);j++){
			delta_y = dy;

			if((j >= (Fy-y_sub_length))
			    &&(j <= (Fy+y_sub_length))){
				delta_y = dy_sub;
			}

			for(k=0;k<=(zs-1);k++){
				hy = hyi[i][j][k]-hyi[i-1][j][k];
				if(i == (Fx-x_sub_length)){
					hy = hyi[i][j][k]-hy_x_absorb[0][j][k];
				}
				if(i == (Fx+x_sub_length)){
					hy = hy_x_absorb[1][j][k]-hyi[i-1][j][k];
				}

				hx = hxi[i][j][k]-hxi[i][j-1][k];
				if(j == (Fy-y_sub_length)){
					hx = hxi[i][j][k]-hx_y_absorb[i][0][k];
				}
				if(j == (Fy+y_sub_length)){
					hx = hx_y_absorb[i][1][k]-hxi[i][j-1][k];
				}