fdtd.cpp

液晶中TE波入射的FDTD算法

					Ex_nxt[i][k] = Ex_pre[i][k]*exp(-1*deltaz[k-398]*dt/EPSILON0)-(1-exp(-1*deltaz[k-398]*dt/EPSILON0))/dz/deltaz[k-398]*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
					Ez_nxt[i][k] = Ez_pre[i][k]+0.5*Z0*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
				}
			}
			//左边
			for (i=1;i<8;i++)
			{
				for (k=8;k<398;k++)
				{
					Hyx_nxt[i][k] = Hyx_pre[i][k]*exp(-1*deltamx[i]*dt/MU0)+(1-exp(-1*deltamx[i]*dt/MU0))/dx/deltamx[i]*(Ez_pre[i+1][k]-Ez_pre[i][k]);
					Hyz_nxt[i][k] = Hyz_pre[i][k]-0.5/Z0*(Ex_pre[i][k+1]-Ex_pre[i][k]);
					Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
					Ex_nxt[i][k] = Ex_pre[i][k]-0.5*Z0*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
					Ez_nxt[i][k] = Ez_pre[i][k]*exp(-1*deltax[i]*dt/EPSILON0)+(1-exp(-1*deltax[i]*dt/EPSILON0))/dx/deltax[i]*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
				}
			}
			//右边
			for (i=1949;i<1956;i++)
			{
				for (k=8;k<398;k++)
				{
					Hyx_nxt[i][k] = Hyx_pre[i][k]*exp(-1*deltamx[i-1948]*dt/MU0)+(1-exp(-1*deltamx[i-1948]*dt/MU0))/dx/deltamx[i-1948]*(Ez_pre[i+1][k]-Ez_pre[i][k]);
					Hyz_nxt[i][k] = Hyz_pre[i][k]-0.5/Z0*(Ex_pre[i][k+1]-Ex_pre[i][k]);
					Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
					Ex_nxt[i][k] = Ex_pre[i][k]-0.5*Z0*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
					Ez_nxt[i][k] = Ez_pre[i][k]*exp(-1*deltax[i-1948]*dt/EPSILON0)+(1-exp(-1*deltax[i-1948]*dt/EPSILON0))/dx/deltax[i-1948]*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
				}
			}

			/************************************************************************/
			/* PML角点(1-7)                                                         */
			/************************************************************************/
			//左下边
			for (i=1;i<8;i++)
			{
				for (k=1;k<8;k++)
				{
					Hyx_nxt[i][k] = Hyx_pre[i][k]*exp(-1*deltamx[i]*dt/MU0)+(1-exp(-1*deltamx[i]*dt/MU0))/dx/deltamx[i]*(Ez_pre[i+1][k]-Ez_pre[i][k]);
					Hyz_nxt[i][k] = Hyz_pre[i][k]*exp(-1*deltamz[k]*dt/MU0)-(1-exp(-1*deltamz[k]*dt/MU0))/dz/deltamz[k]*(Ex_pre[i][k+1]-Ex_pre[i][k]);
					Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
					Ex_nxt[i][k] = Ex_pre[i][k]*exp(-1*deltaz[k]*dt/EPSILON0)-(1-exp(-1*deltaz[k]*dt/EPSILON0))/dz/deltaz[k]*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
					Ez_nxt[i][k] = Ez_pre[i][k]*exp(-1*deltax[i]*dt/EPSILON0)+(1-exp(-1*deltax[i]*dt/EPSILON0))/dx/deltax[i]*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
				}
			}
			//左上边
			for (i=1;i<8;i++)
			{
				for (k=398;k<405;k++)
				{
					Hyx_nxt[i][k] = Hyx_pre[i][k]*exp(-1*deltamx[i]*dt/MU0)+(1-exp(-1*deltamx[i]*dt/MU0))/dx/deltamx[i]*(Ez_pre[i+1][k]-Ez_pre[i][k]);
					Hyz_nxt[i][k] = Hyz_pre[i][k]*exp(-1*deltamz[k-397]*dt/MU0)-(1-exp(-1*deltamz[k-397]*dt/MU0))/dz/deltamz[k-397]*(Ex_pre[i][k+1]-Ex_pre[i][k]);
					Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
					Ex_nxt[i][k] = Ex_pre[i][k]*exp(-1*deltaz[k-397]*dt/EPSILON0)-(1-exp(-1*deltaz[k-397]*dt/EPSILON0))/dz/deltaz[k-397]*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
					Ez_nxt[i][k] = Ez_pre[i][k]*exp(-1*deltax[i]*dt/EPSILON0)+(1-exp(-1*deltax[i]*dt/EPSILON0))/dx/deltax[i]*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
				}
			}
			//右上边
			for (i=1948;i<1955;i++)
			{
				for (k=398;k<405;k++)
				{
					Hyx_nxt[i][k] = Hyx_pre[i][k]*exp(-1*deltamx[i-1947]*dt/MU0)+(1-exp(-1*deltamx[i-1947]*dt/MU0))/dx/deltamx[i-1937]*(Ez_pre[i+1][k]-Ez_pre[i][k]);
					Hyz_nxt[i][k] = Hyz_pre[i][k]*exp(-1*deltamz[k-397]*dt/MU0)-(1-exp(-1*deltamz[k-397]*dt/MU0))/dz/deltamz[k-397]*(Ex_pre[i][k+1]-Ex_pre[i][k]);
					Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
					Ex_nxt[i][k] = Ex_pre[i][k]*exp(-1*deltaz[k-397]*dt/EPSILON0)-(1-exp(-1*deltaz[k-397]*dt/EPSILON0))/dz/deltaz[k-397]*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
					Ez_nxt[i][k] = Ez_pre[i][k]*exp(-1*deltax[i-1947]*dt/EPSILON0)+(1-exp(-1*deltax[i-1937]*dt/EPSILON0))/dx/deltax[i-1937]*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
				}
			}
			//右下边
			for (i=1948;i<1955;i++)
			{
				for (k=1;k<8;k++)
				{
					Hyx_nxt[i][k] = Hyx_pre[i][k]*exp(-1*deltamx[i-1947]*dt/MU0)+(1-exp(-1*deltamx[i-1947]*dt/MU0))/dx/deltamx[i-1947]*(Ez_pre[i+1][k]-Ez_pre[i][k]);
					Hyz_nxt[i][k] = Hyz_pre[i][k]*exp(-1*deltamz[k]*dt/MU0)-(1-exp(-1*deltamz[k]*dt/MU0))/dz/deltamz[k]*(Ex_pre[i][k+1]-Ex_pre[i][k]);
					Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
					Ex_nxt[i][k] = Ex_pre[i][k]*exp(-1*deltaz[k]*dt/EPSILON0)-(1-exp(-1*deltaz[k]*dt/EPSILON0))/dz/deltaz[k]*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
					Ez_nxt[i][k] = Ez_pre[i][k]*exp(-1*deltax[i-1947]*dt/EPSILON0)+(1-exp(-1*deltax[i-1947]*dt/EPSILON0))/dx/deltax[i-1947]*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
				}
			}

			/************************************************************************/
			/* PML外边界点(0)                                                       */
			/************************************************************************/
			//下边
			for (i=0;i<1956;i++)
			{
				Hyx_nxt[i][0] = Hyx_pre[i][0]+0.5/Z0*(Ez_pre[i+1][0]-Ez_pre[i][0]);
				Hyz_nxt[i][0] = Hyz_pre[i][0]*exp(-1*deltamz[0]*dt/MU0)-(1-exp(-1*deltamz[0]*dt/MU0))/dz/deltamz[0]*(Ex_pre[i][1]-Ex_pre[i][0]);
				Hy_nxt[i][0] = Hyz_nxt[i][0]+Hyz_nxt[i][0];
				Ex_nxt[i][0] = Ex_pre[i][0]*exp(-1*deltaz[0]*dt/EPSILON0)-(1-exp(-1*deltaz[0]*dt/EPSILON0))/dz/deltaz[0]*Hy_nxt[i][0];
				Ez_nxt[i][0] = Ez_pre[i][0]+0.5*Z0*Hy_nxt[i][0];
			//上边
				Hyx_nxt[i][405] = Hyx_pre[i][405]+0.5/Z0*(0-Ez_pre[i][405]);
				Hyz_nxt[i][405] = Hyz_pre[i][405]*exp(-1*deltamz[0]*dt/MU0)-(1-exp(-1*deltamz[0]*dt/MU0))/dz/deltamz[0]*(0-Ex_pre[i][405]);
				Hy_nxt[i][405] = Hyz_nxt[i][405]+Hyz_nxt[i][405];
				Ex_nxt[i][405] = Ex_pre[i][405]*exp(-1*deltaz[0]*dt/EPSILON0)-(1-exp(-1*deltaz[0]*dt/EPSILON0))/dz/deltaz[0]*(Hy_nxt[i][405]-Hy_nxt[i][404]);
				Ez_nxt[i][405] = Ez_pre[i][405]+0.5*Z0*(Hy_nxt[i][405]-Hy_nxt[i-1][405]);
			}
			//左边
			for (k=0;k<406;k++)
			{
					Hyx_nxt[0][k] = Hyx_pre[0][k]*exp(-1*deltamx[0]*dt/MU0)+(1-exp(-1*deltamx[0]*dt/MU0))/dx/deltamx[0]*(Ez_pre[1][k]-Ez_pre[0][k]);
					Hyz_nxt[0][k] = Hyz_pre[0][k]-0.5/Z0*(Ex_pre[0][k+1]-Ex_pre[0][k]);
					Hy_nxt[0][k] = Hyz_nxt[0][k]+Hyz_nxt[0][k];
					Ex_nxt[0][k] = Ex_pre[0][k]-0.5*Z0*(Hy_nxt[0][k]-Hy_nxt[0][k-1]);
					Ez_nxt[0][k] = Ez_pre[0][k]*exp(-1*deltax[0]*dt/EPSILON0)+(1-exp(-1*deltax[0]*dt/EPSILON0))/dx/deltax[0]*Hy_nxt[0][k];
			//右边
					Hyx_nxt[1955][k] = Hyx_pre[1955][k]*exp(-1*deltamx[0]*dt/MU0)+(1-exp(-1*deltamx[0]*dt/MU0))/dx/deltamx[0]*(0-Ez_pre[1955][k]);
					Hyz_nxt[1955][k] = Hyz_pre[1955][k]-0.5/Z0*(Ex_pre[1955][k+1]-Ex_pre[1955][k]);
					Hy_nxt[1955][k] = Hyz_nxt[1955][k]+Hyz_nxt[1955][k];
					Ex_nxt[1955][k] = Ex_pre[1955][k]-0.5*Z0*(Hy_nxt[1955][k]-Hy_nxt[1955][k-1]);
					Ez_nxt[1955][k] = Ez_pre[1955][k]*exp(-1*deltax[0]*dt/EPSILON0)+(1-exp(-1*deltax[0]*dt/EPSILON0))/dx/deltax[0]*(Hy_nxt[1955][k]-Hy_nxt[1954][k]);
			}
			/************************************************************************/
			/*真空中FDTD公式                                                        */
			/************************************************************************/
			for (i=8;i<1948;i++)   
			{
				for (k=8;k<398;k++)
				{
					Hy_nxt[i][k] = Hy_pre[i][k]+dt*K0*((Ez_pre[i+1][k]-Ez_pre[i][k])/dx-(Ex_pre[i][k+1]-Ex_pre[i][k])/dz);
					Ex_nxt[i][k] = Ex_pre[i][k]-dt*K0*(Hy_nxt[i][k]-Hy_nxt[i][k-1])/dz;	
					Ez_nxt[i][k] = Ez_pre[i][k]+dt*K0*(Hy_nxt[i][k]-Hy_nxt[i-1][k])/dx;
				}
			}
			
			/************************************************************************/
			/*液晶层中FDTD公式                                                      */
			/************************************************************************/
			for (i=78;i<1878;i++) 
			{
				for (k=78;k<328;k++)
				{
					temp1 = 13.8*sin(thita[i-78][k-78])*cos(thita[i-78][k-78]);
					temp2 = 19.0+13.8*pow(cos(thita[i-78][k-78]),2);
					temp3 = 19.0+13.8*pow(sin(thita[i-78][k-78]),2);
					temp4 = pow(temp1,2)-temp2*temp3;
					Ex_nxt[i][k] = Ex_pre[i][k]+dt*temp1/temp4*K0*(Hy_nxt[i][k]-Hy_nxt[i-1][k])/dx+dt*temp3/temp4*K0*(Hy_nxt[i][k]-Hy_nxt[i][k-1])/dz;
					Ez_nxt[i][k] = Ez_pre[i][k]-dt*temp2/temp4*K0*(Hy_nxt[i][k]-Hy_nxt[i-1][k])/dx-dt*temp1/temp4*K0*(Hy_nxt[i][k]-Hy_nxt[i][k-1])/dz;
				}
			}

			/************************************************************************/
			/*真空与液晶连接边界                                                    */
			/************************************************************************/
			for (i=78;i<1878;i++)
			{
				//下边界
				Ex_nxt[i][78] = (-1)*Ex_nxt[i][79]+EPSILON0/epsilonzxf[i-78]*(Ez_nxt[i][77]+Ez_nxt[i+1][77])-epsilonzzf[i-78]/epsilonzxf[i-78]*(Ez_nxt[i][78]+Ez_nxt[i+1][78]);
				//上边界
				Ex_nxt[i][328] = (-1)*Ex_nxt[i][327]+EPSILON0/epsilonzxb[i-78]*(Ez_nxt[i][328]+Ez_nxt[i+1][328])-epsilonzzb[i-78]/epsilonzxb[i-78]*(Ez_nxt[i][327]+Ez_nxt[i+1][327]);
			}
			for (k=78;k<328;k++)
			{
				//左边界
				Ez_nxt[78][k] = (-1)*Ez_nxt[78][k]+EPSILON0/epsilonxzl[k-78]*(Ex_nxt[77][k]+Ex_nxt[77][k+1])-epsilonxxl[k-78]/epsilonxzl[k-78]*(Ex_nxt[78][k]+Ex_nxt[78][k+1]);
				//右边界
				Ez_nxt[1878][k] = (-1)*Ez_nxt[1877][k]+EPSILON0/epsilonxzr[k-78]*(Ex_nxt[1878][k]+Ex_nxt[1878][k+1])-epsilonxxr[k-78]/epsilonxzr[k-78]*(Ex_nxt[1877][k]+Ex_nxt[1877][k+1]);
			}
			
			/************************************************************************/
			/* 总场边界 (包含入射波的引入)                                          */
			/************************************************************************/
			//下边界
			E0 = 7.74;
		    Exincf_nxt = E0*sin(omega*steps*dt-2*PI/lamda*(58-m0)*dz);
			Hyincf_nxt = E0/Z0*sin(omega*(steps-0.5)*dt-2*PI/lamda*(58.5-m0)*dz);
			Exincf_pre = E0*sin(omega*(steps-1)*dt-2*PI/lamda*(58-m0)*dz);
			Hyincf_pre = E0/Z0*sin(omega*(steps-1.5)*dt-2*PI/lamda*(58.5-m0)*dz);
			//for (i=59;i<1897;i++)
			for (i=8;i<1948;i++)
			{
				Hy_nxt[i][58] = Hy_pre[i][58]+dt*K0*((Ez_pre[i+1][58]-Ez_pre[i][58])/dx-(Ex_pre[i][59]-Ex_pre[i][58])/dz)-0.5*Hyincf_nxt;
				Ex_nxt[i][58] = Ex_pre[i][58]-dt*K0/dz*(Hy_nxt[i][58]-Hy_nxt[i][57])-0.5*Exincf_nxt;
			}
			Hy_nxt[58][58] = Hy_pre[58][58]+dt*K0*((Ez_pre[59][58]-Ez_pre[58][58])/dx-(Ex_pre[58][59]-Ex_pre[58][58])/dz)-0.5*Hyincf_nxt;
			Ex_nxt[58][58] = Ex_pre[58][58]-dt*K0/dz*(Hy_nxt[58][58]-Hy_nxt[58][57])-0.5*Exincf_nxt;
			Hy_nxt[1897][58] = Hy_pre[1897][58]+dt*K0*((Ez_pre[1898][58]-Ez_pre[1897][58])/dx-(Ex_pre[1897][59]-Ex_pre[1897][58])/dz)-0.5*Hyincf_nxt;
			Ex_nxt[1897][58] = Ex_pre[1897][58]-dt*K0/dz*(Hy_nxt[1897][58]-Hy_nxt[1897][57])-0.5*Exincf_nxt;
			/*上边界
			//Hyincb = sqrt(EPSILON0/MU0)*E0*sin(omega*steps*dt-2*PI/lamda*(348-m0)*dz);
			Hyincb = 0;
			//Exincb = E0*sin(omega*(steps-1)*dt-2*PI/lamda*(348-m0)*dz);
			Exincb = 0;
			for (i=58;i<1898;i++)
			{
				Hy_nxt[i][348] = Hy_pre[i][348]+dt/Z0*((Ez_pre[i+1][348]-Ez_pre[i][348])/dx-(Ex_pre[i][349]-Ex_pre[i][348]+Exincb)/dz);
				Ex_nxt[i][348] = Ex_pre[i][348]-dt*Z0/dz*(Hy_nxt[i][348]-Hy_nxt[i][347]+Hyincb);
			}
			//左边界
			for (k=58;k<348;k++)
			{
				//Hyinc[k-58] = sqrt(EPSILON0/MU0)*E0*sin(omega*steps*dt-2*PI/lamda*(k-m0)*dz);
				Hyinc[k-58] = 0;
				//Ezinc[k-58] = 0;
				Hy_nxt[57][k] = Hy_pre[57][k]+dt/Z0*((Ez_pre[58][k]-Ez_pre[57][k])/dx-(Ex_pre[57][k+1]-Ex_pre[57][k])/dz);
				Ez_nxt[58][k] = Ez_pre[58][k]+dt*Z0/dx*(Hy_nxt[58][k]-Hy_nxt[57][k]-Hyinc[k-58]);
			}
			//右边界
			for (k=58;k<348;k++)
			{
				Hy_nxt[1898][k] = Hy_pre[1898][k]+dt/Z0*((Ez_pre[1899][k]-Ez_pre[1898][k])/dx-(Ex_pre[1898][k+1]-Ex_pre[1898][k])/dz);
				Ez_nxt[1898][k] = Ez_pre[1898][k]+dt*Z0/dx*(Hy_nxt[1898][k]-Hy_nxt[1897][k]+Hyinc[k-58]);
			}*/			
			
			/************************************************************************/                                                    
			/* 矩阵值传递                                                           */
			/************************************************************************/
			for (i=0;i<NUM_EX;i++)
			{
				for (k=0;k<NUM_EZ;k++)
				{
					Ex_pre[i][k] = Ex_nxt[i][k];
					Ez_pre[i][k] = Ez_nxt[i][k];
					Hy_pre[i][k] = Hy_nxt[i][k];
				}
			}
			ShowProInfo();
		}

		/************************************************************************/
		/* FDTD计算完毕                                                         */
		/************************************************************************/
		for (i=0;i<NUM_EX;i++)
			{
				for (k=0;k<NUM_EZ;k++)
				{
					Ex_final[i+k*NUM_EX] = Ex_nxt[i][k];
					Ez_final[i+k*NUM_EX] = Ez_nxt[i][k];
					Hy_final[i+k*NUM_EX] = Hy_nxt[i][k];
				}
			}
	//}	
	if (WriteMatFile()==1)
	{
		return(EXIT_FAILURE);
	}

	printf("Done!\n");
	return(EXIT_SUCCESS);
}