part_4.br

来自「The Finite Difference Time Domain Method」· BR 代码 · 共 83 行

// Initialize PML Boundaries
for (i=1; i<1+iPML; i++) {

	ii=(float) i;
	ii=(1+iPML)-ii;
		
	sig1=pow(((ii-0.5)/iPML),m)*sigmax;
	sig2=(mu0/eps0)*pow(((ii)/iPML),m)*sigmax;

	a1=pow(((iPML-(ii-1+.5))/iPML),m)*amax;
	a2=(mu0/eps0)*pow(((iPML-(ii-1))/iPML),m)*amax;
			
	k1=1+(kappa-1)*pow(((ii-0.5)/iPML),m);
	k2=1+(kappa-1)*pow(((ii)/iPML),m);

	aKex[i+1]=1/k1;
	aKhx[i]=1/k2;	

	aKex[nx-i-1]=1/k1;
	aKhx[nx-i-1]=1/k2;	

	abex[i+1]=exp((-dt/eps0)*((sig1/k1)+a1));
	acex[i+1]=((sig1/dx)/((sig1*k1)+(a1*k1*k1))) *
		  (exp((-dt/eps0)*((sig1/k1)+a1))-1);


	abex[nx-i-1]=exp((-dt/eps0)*((sig1/k1)+a1));
	acex[nx-i-1]=((sig1/dx)/((sig1*k1)+(a1*k1*k1))) *
		     (exp((-dt/eps0)*((sig1/k1)+a1))-1);

	abhx[i]=exp((-dt/mu0)*((sig2/k2)+a2));
	achx[i]=((sig2/dx)/((sig2*k2)+(a2*k2*k2))) *
	        (exp((-dt/mu0)*((sig2/k2)+a2))-1);	

	abhx[nx-i-1]=exp((-dt/mu0)*((sig2/k2)+a2));
	achx[nx-i-1]=((sig2/dx)/((sig2*k2)+(a2*k2*k2))) *
		     (exp((-dt/mu0)*((sig2/k2)+a2))-1);

	aKey[i+1]=1/k1;
	aKhy[i]=1/k2;	

	aKey[ny-i-1]=1/k1;
	aKhy[ny-i-1]=1/k2;	

	abey[i+1]=exp((-dt/eps0)*((sig1/k1)+a1));
	acey[i+1]=((sig1/dx)/((sig1*k1)+(a1*k1*k1))) *
		  (exp((-dt/eps0)*((sig1/k1)+a1))-1);


	abey[ny-i-1]=exp((-dt/eps0)*((sig1/k1)+a1));
	acey[ny-i-1]=((sig1/dx)/((sig1*k1)+(a1*k1*k1))) *
		     (exp((-dt/eps0)*((sig1/k1)+a1))-1);

	abhy[i]=exp((-dt/mu0)*((sig2/k2)+a2));
	achy[i]=((sig2/dx)/((sig2*k2)+(a2*k2*k2))) *
		(exp((-dt/mu0)*((sig2/k2)+a2))-1);	

	abhy[ny-i-1]=exp((-dt/mu0)*((sig2/k2)+a2));
	achy[ny-i-1]=((sig2/dx)/((sig2*k2)+(a2*k2*k2))) *
		     (exp((-dt/mu0)*((sig2/k2)+a2))-1);
}

	
	
	// Place Our Point Source (Placed in the approximate center)
 	aCs[nx*(ysize/2)+(xsize/2)]=1/dx;

	// Place Sample Scattering Object Er=10.2

        // This places a rectangular scattering 
        // Dielectric box 10 cells from the point source
        // That is 16 cells wide and 10 Cells Deep

	for (j=((ysize/2)-8); j<((ysize/2)+9); j++) {
		for (i=((xsize/2)+10); i<((xsize/2)+20); i++) {
			aCeze[nx*j+i]=1;
			aCezh[nx*j+i]=(2*ce)/(10.2);
		}
	}

	fprintf (pFile,"Cezh=%e Chxe=%e\n",aCezh[10],aChxe[10]);