📄 nf.cpp
字号:
//PML
#include "my.h"
int main()
{
Initial();
int timeloop;
int i,j,k;
for (timeloop=0;timeloop<=NN;timeloop++)
{
for(i=1;i<=Nx-1;i++)
{
for(j=0;j<=Ny-1;j++)
{
for(k=0;k<=Nz-1;k++)
{
if(i==mx&&(k==mz-1||k==mz)&&j<=myur&&j>=mydl)
{
Hxz[i][j][k]=Hxz[i][j][k]
+dt/miu*2/log(dy/r)*(Ey[i][j][k+1]-Ey[i][j][k])/dz;
}
else
{
Hxy[i][j][k]=(1-sgmmy[j]*dt/(2*miu))/(1+sgmmy[j]*dt/(2*miu))*Hxy[i][j][k]+(dt/miu)/(1+sgmmy[j]*dt/(2*miu))*(Ez[i][j][k]-Ez[i][j+1][k])/dy;
Hxz[i][j][k]=(1-sgmmz[k]*dt/(2*miu))/(1+sgmmz[k]*dt/(2*miu))*Hxz[i][j][k]+(dt/miu)/(1+sgmmz[k]*dt/(2*miu))*(Ey[i][j][k+1]-Ey[i][j][k])/dz;
Hx[i][j][k]=Hxy[i][j][k]+Hxz[i][j][k];
}
}
}
}
for(i=0;i<=Nx-1;i++)
{
for(j=1;j<=Ny-1;j++)
{
for( k=0;k<=Nz-1;k++)
{
if(i==mx-1&&j>=mydl&&j<=myur&&k==mz-1)
{
Hyx[i][j][k]=Hyx[i][j][k]+dt/miu/(dx*dz-pi*r*r/4.0)*dz
*(log(dz/r)/2.0*Ez[i+1][j][k]-Ez[i][j][k]);
Hyz[i][j][k]=Hyz[i][j][k]-dt/miu/(dx*dz-pi*r*r/4.0)*dx
*(log(dx/r)/2.0*Ex[i][j][k+1]-Ex[i][j][k]);
}
else if(i==mx&&j<=myur&&j>=mydl&&k==mz)
{
Hyx[i][j][k]=Hyx[i][j][k]+dt/miu/(dx*dz-pi*r*r/4.0)*dz
*(Ez[i+1][j][k]-log(dz/r)/2.0*Ez[i][j][k]);
Hyz[i][j][k]=Hyz[i][j][k]-dt/miu/(dx*dz-pi*r*r/4.0)*dx
*(Ex[i][j][k+1]-log(dx/r)/2.0*Ex[i][j][k]);
}
else if(i==mx&&j<=myur&&j>=mydl&&k==mz-1)
{
Hyx[i][j][k]=Hyx[i][j][k]+dt/miu/(dx*dz-pi*r*r/4.0)*dz
*(Ez[i+1][j][k]-log(dz/r)/2.0*Ez[i][j][k]);
Hyz[i][j][k]=Hyz[i][j][k]-dt/miu/(dx*dz-pi*r*r/4.0)*dx
*(log(dx/r)/2.0*Ex[i][j][k+1]-Ex[i][j][k]);
}
else if(i==mx-1&&j<=myur&&j>=mydl&&k==mz)
{
Hyx[i][j][k]=Hyx[i][j][k]+dt/miu/(dx*dz-pi*r*r/4.0)*dz
*(log(dz/r)/2.0*Ez[i+1][j][k]-Ez[i][j][k]);
Hyz[i][j][k]=Hyz[i][j][k]-dt/miu/(dx*dz-pi*r*r/4.0)*dx
*(Ex[i][j][k+1]-log(dx/r)/2.0*Ex[i][j][k]);
}
else
{
Hyx[i][j][k]=(1-sgmmx[i]*dt/(2*miu))/(1+sgmmx[i]*dt/(2*miu))*Hyx[i][j][k]+(dt/miu)/(1+sgmmx[i]*dt/(2*miu))*(Ez[i+1][j][k]-Ez[i][j][k])/dx;
Hyz[i][j][k]=(1-sgmmz[k]*dt/(2*miu))/(1+sgmmz[k]*dt/(2*miu))*Hyz[i][j][k]+(dt/miu)/(1+sgmmz[k]*dt/(2*miu))*(Ex[i][j][k]-Ex[i][j][k+1])/dz;
}
}
}
}
for( i=0;i<=Nx-1;i++)
{
for( j=0;j<=Ny-1;j++)
{
for( k=0;k<=Nz-1;k++)
{
if((i==mx||i==mx-1)&&k==mz&&j<=myur&&j>=mydl)
{
Hzx[i][j][k]=Hzx[i][j][k]
-dt/miu*2.0/log(dx/r)*(Ey[i+1][j][k]-Ey[i][j][k])/dx;
}
else
{
Hzx[i][j][k]=(1-sgmmx[i]*dt/(2*miu))/(1+sgmmx[i]*dt/(2*miu))*Hzx[i][j][k]+(dt/miu)/(1+sgmmx[i]*dt/(2*miu))*(Ey[i][j][k]-Ey[i+1][j][k])/dx;
Hzy[i][j][k]=(1-sgmmy[j]*dt/(2*miu))/(1+sgmmy[j]*dt/(2*miu))*Hzy[i][j][k]+(dt/miu)/(1+sgmmy[j]*dt/(2*miu))*(Ex[i][j+1][k]-Ex[i][j][k])/dy;
}
Hz[i][j][k]=Hzx[i][j][k]+Hzy[i][j][k];
}
}
}
for(i=0;i<=Nx-1;i++)
{
for(j=1;j<=Ny-1;j++)
{
for(k=1;k<=Nz-1;k++)
{
Exy[i][j][k]=(1-sgmy[j]*dt/(2*ee))/(1+sgmy[j]*dt/(2*ee))*Exy[i][j][k]+(dt/ee)/(1+sgmy[j]*dt/(2*ee))*(Hzx[i][j][k]+Hzy[i][j][k]-Hzx[i][j-1][k]-Hzy[i][j-1][k])/dy;
Exz[i][j][k]=(1-sgmz[k]*dt/(2*ee))/(1+sgmz[k]*dt/(2*ee))*Exz[i][j][k]+(dt/ee)/(1+sgmz[k]*dt/(2*ee))*(Hyz[i][j][k-1]+Hyx[i][j][k-1]-Hyz[i][j][k]-Hyx[i][j][k])/dz;
Ex[i][j][k]=Exy[i][j][k]+Exz[i][j][k];
}
}
}
for(i=1;i<=Nx-1;i++)
{
for(j=0;j<=Ny-1;j++)
{
for(k=1;k<=Nz-1;k++)
{
if(i==mx&&k==mz&&j<=myur&&j>=mydl)
{
Eyx[i][j][k]=0;
Eyz[i][j][k]=0;
}
else
{
Eyx[i][j][k]=(1-sgmx[i]*dt/(2*ee))/(1+sgmx[i]*dt/(2*ee))*Eyx[i][j][k]+(dt/ee)/(1+sgmx[i]*dt/(2*ee))*(Hzx[i-1][j][k]+Hzy[i-1][j][k]-Hzx[i][j][k]-Hzy[i][j][k])/dx;
Eyz[i][j][k]=(1-sgmz[k]*dt/(2*ee))/(1+sgmz[k]*dt/(2*ee))*Eyz[i][j][k]+(dt/ee)/(1+sgmz[k]*dt/(2*ee))*(Hxy[i][j][k]+Hxz[i][j][k]-Hxy[i][j][k-1]-Hxz[i][j][k-1])/dz;
}
Ey[i][j][k]=Eyx[i][j][k]+Eyz[i][j][k];
}
}
}
for(i=1;i<=Nx-1;i++)
{
for(j=1;j<=Ny-1;j++)
{
for(k=0;k<=Nz-1;k++)
{
Ezx[i][j][k]=(1-sgmx[i]*dt/(2*ee))/(1+sgmx[i]*dt/(2*ee))*Ezx[i][j][k]+(dt/ee)/(1+sgmx[i]*dt/(2*ee))*(Hyz[i][j][k]+Hyx[i][j][k]-Hyz[i-1][j][k]-Hyx[i-1][j][k])/dx;
Ezy[i][j][k]=(1-sgmy[j]*dt/(2*ee))/(1+sgmy[j]*dt/(2*ee))*Ezy[i][j][k]+(dt/ee)/(1+sgmy[j]*dt/(2*ee))*(Hxy[i][j-1][k]+Hxz[i][j-1][k]-Hxy[i][j][k]-Hxz[i][j][k])/dy;
Ez[i][j][k]=Ezx[i][j][k]+Ezy[i][j][k];
}
}
}
Ey[Nx/2][Ny/2][Nz/2]=Ey[Nx/2][Ny/2][Nz/2]+sin(2*pi*f*((timeloop+0.5)*dt))/Gap;
E1[timeloop]=Exy[Nx/2][Ny/2][Nz/2+8]+Exz[Nx/2][Ny/2][Nz/2+8];
E2[timeloop]=Exy[Nx/2][Ny/2][Nz/2-8]+Exz[Nx/2][Ny/2][Nz/2-8];
printf(" timestep:%d \n",timeloop);
if(timeloop>=NN-120)
cal_J(timeloop);
}
fangxiangtu();
return 1;
}
void Initial()
{
int loop1,loop2,loop3;
for( int loop=0;loop<=Nx-1;loop++)
{
if(loop<=NNN-1)
sgmmx[loop]=sgmmmax*(NNN-1.0/2.0-loop)/NNN*(NNN-1.0/2.0-loop)/NNN;
else if(loop>=Nx-NNN)
sgmmx[loop]=sgmmmax*(loop-Nx+NNN+1.0/2.0)/NNN*(loop-Nx+NNN+1.0/2.0)/NNN;
else
sgmmx[loop]=0;
}
for(loop=0;loop<=Ny-1;loop++)
{
if(loop<=NNN-1)
sgmmy[loop]=sgmmmay*(NNN-1.0/2.0-loop)/NNN*(NNN-1.0/2.0-loop)/NNN;
else if(loop>=Ny-NNN)
sgmmy[loop]=sgmmmay*(loop-Ny+NNN+1.0/2.0)/NNN*(loop-Ny+NNN+1.0/2.0)/NNN;
else
sgmmy[loop]=0;
}
for(loop=0;loop<=Nz-1;loop++)
{
if(loop<=NNN-1)
sgmmz[loop]=sgmmmaz*(NNN-1.0/2.0-loop)/NNN*(NNN-1.0/2.0-loop)/NNN;
else if(loop>=Nz-NNN)
sgmmz[loop]=sgmmmaz*(loop-Nz+NNN+1.0/2.0)/NNN*(loop-Nz+NNN+1.0/2.0)/NNN;
else
sgmmz[loop]=0;
}
for(loop=0;loop<=Nx;loop++)
{
if(loop<=NNN)
sgmx[loop]=sgmmax*(NNN-loop)/NNN*(NNN-loop)/NNN;
else if(loop>=Nx-NNN)
sgmx[loop]=sgmmax*(loop-Nx+NNN)/NNN*(loop-Nx+NNN)/NNN;
else
sgmx[loop]=0;
}
for(loop=0;loop<=Ny;loop++)
{
if(loop<=NNN)
sgmy[loop]=sgmmay*(NNN-loop)/NNN*(NNN-loop)/NNN;
else if(loop>=Ny-NNN)
sgmy[loop]=sgmmay*(loop-Ny+NNN)/NNN*(loop-Ny+NNN)/NNN;
else
sgmy[loop]=0;
}
for(loop=0;loop<=Nz;loop++)
{
if(loop<=NNN)
sgmz[loop]=sgmmaz*(NNN-loop)/NNN*(NNN-loop)/NNN;
else if(loop>=Nz-NNN)
sgmz[loop]=sgmmaz*(loop-Nz+NNN)/NNN*(loop-Nz+NNN)/NNN;
else
sgmz[loop]=0;
}
for(loop1=0;loop1<Nx+1;loop1++)
for(loop2=0;loop2<Ny+1;loop2++)
for(loop3=0;loop3<Nz;loop3++)
{ Ezx[loop1][loop2][loop3]=0;Ezy[loop1][loop2][loop3]=0;}
for(loop1=0;loop1<Nx;loop1++)
for(loop2=0;loop2<Ny+1;loop2++)
for(loop3=0;loop3<Nz+1;loop3++)
{ Exy[loop1][loop2][loop3]=0; Exz[loop1][loop2][loop3]=0;}
for(loop1=0;loop1<Nx+1;loop1++)
for(loop2=0;loop2<Ny;loop2++)
for(loop3=0;loop3<Nz+1;loop3++)
{ Eyx[loop1][loop2][loop3]=0; Eyz[loop1][loop2][loop3]=0;}
for(loop1=0;loop1<Nx;loop1++)
for(loop2=0;loop2<Ny;loop2++)
for(loop3=0;loop3<Nz+1;loop3++)
{ Hzx[loop1][loop2][loop3]=0; Hzy[loop1][loop2][loop3]=0;}
for(loop1=0;loop1<Nx+1;loop1++)
for(loop2=0;loop2<Ny;loop2++)
for(loop3=0;loop3<Nz;loop3++)
{ Hxy[loop1][loop2][loop3]=0; Hxz[loop1][loop2][loop3]=0;}
for(loop1=0;loop1<Nx;loop1++)
for(loop2=0;loop2<Ny+1;loop2++)
for(loop3=0;loop3<Nz;loop3++)
{ Hyx[loop1][loop2][loop3]=0; Hyz[loop1][loop2][loop3]=0;}
for(loop1=0;loop1<NN/2;loop1++)
{
pp1[loop1].shi=0;pp1[loop1].xu=0;
pp2[loop1].shi=0;pp2[loop1].xu=0;
zk[loop1].shi=0;zk[loop1].xu=0;
}
fxhs1.sita.shi=0;fxhs1.sita.xu=0;
fxhs1.fai.shi=0;fxhs1.fai.xu=0;
for(loop2=0;loop2<180;loop2++)
{fxhs[loop2].shi=0;fxhs[loop2].xu=0;}
for(loop1=0;loop1<180;loop1++)
{si[loop1]=sin(loop1/180.0*pi);co[loop1]=cos(loop1/180.0*pi);}
for(loop1=0;loop1<Nx-2*NNN-2*DDD;loop1++)
for(loop2=0;loop2<Ny-2*NNN-2*DDD+1;loop2++)
{
Jxu[loop1][loop2].zhengfu=0; Jxu[loop1][loop2].xiangwei=0;
Jxd[loop1][loop2].zhengfu=0; Jxd[loop1][loop2].xiangwei=0;
Jmxu[loop1][loop2].zhengfu=0; Jmxu[loop1][loop2].xiangwei=0;
Jmxd[loop1][loop2].zhengfu=0; Jmxd[loop1][loop2].xiangwei=0;
}
for(loop1=0;loop1<Nx-2*NNN-2*DDD;loop1++)
for(loop2=0;loop2<Nz-2*NNN-2*DDD+1;loop2++)
{
Jmxl[loop1][loop2].zhengfu=0; Jmxl[loop1][loop2].xiangwei=0;
Jmxr[loop1][loop2].zhengfu=0; Jmxr[loop1][loop2].xiangwei=0;
Jxl[loop1][loop2].zhengfu=0; Jxl[loop1][loop2].xiangwei=0;
Jxr[loop1][loop2].zhengfu=0; Jxr[loop1][loop2].xiangwei=0;
}
for(loop1=0;loop1<Nx-2*NNN-2*DDD;loop1++)
for(loop2=0;loop2<Nz-2*NNN-2*DDD+1;loop2++)
{
Jzl[loop1][loop2].zhengfu=0; Jzl[loop1][loop2].xiangwei=0;
Jzr[loop1][loop2].zhengfu=0; Jzr[loop1][loop2].xiangwei=0;
Jmzl[loop1][loop2].zhengfu=0; Jmzl[loop1][loop2].xiangwei=0;
Jmzr[loop1][loop2].zhengfu=0; Jmzr[loop1][loop2].xiangwei=0;
}
for(loop1=0;loop1<Nx-2*NNN-2*DDD+1;loop1++)
for(loop2=0;loop2<Ny-2*NNN-2*DDD;loop2++)
{
Jyu[loop1][loop2].zhengfu=0; Jyu[loop1][loop2].xiangwei=0;
Jyd[loop1][loop2].zhengfu=0; Jyd[loop1][loop2].xiangwei=0;
Jmyu[loop1][loop2].zhengfu=0; Jmyu[loop1][loop2].xiangwei=0;
Jmyd[loop1][loop2].zhengfu=0; Jmyd[loop1][loop2].xiangwei=0;
}
for(loop1=0;loop1<Ny-2*NNN-2*DDD+1;loop1++)
for(loop2=0;loop2<Nz-2*NNN-2*DDD;loop2++)
{
Jyf[loop1][loop2].zhengfu=0; Jyf[loop1][loop2].xiangwei=0;
Jyb[loop1][loop2].zhengfu=0; Jyb[loop1][loop2].xiangwei=0;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -