📄 xiangyang26.cpp
字号:
/*PML完全匹配层*/
#include "shuju7.h"
#include "cbb.h"
#include "bianjie.h"
main()
{ FILE *out;
char outfile[3];
int m,lp,loop,loop1,i,j,k,start=1000;
float Iz1,Uz1,Uzz1[NN],Izz1[NN],Iz2,Uz2,Uzz2[NN],Izz2[NN],fs=0;
printf("%f,%f,%f,%f,%f,%f\n",dx,dy,dz,r0,r1,c/f);
printf("%f\n",1/sqrt((1/dx)*(1/dx)+(1/dy)*(1/dy)+(1/dz)*(1/dz)));
printf("%f\n",1/(dx*dy/dz));
// fuzhi();
// printf("%f\n",1.0/(2*Pi*f*dt));
for(DD1=20;DD1<21;DD1=DD1+2)
for(DD2=20;DD2<21;DD2=DD2+2)
// for (movedis=-10;movedis<=10;movedis=movedis+5)
// {
// fuzhi();
// for(lp=0;lp<10;lp++)
// { r0=(c/f)/400.0;r1=(c/f)/(50+20.0*lp);
// printf("%f\n",1/sqrt((1/dx)*(1/dx)+(1/dy)*(1/dy)+(1/dz)*(1/dz)));
// printf("%f\n",1/(2*Pi*f*ee*er*dx*dy/dz));
{ fuzhi();
printf("dfskldfk");
for(loop=1;loop<=NN;loop++) /*主时间循环*/
{
m=chuanbo1(loop);
m=bianjie();
// Uz=-exp(-(loop*dt-100*dt)*(loop*dt-100*dt)/(dt*dt*1000))*sin(2*Pi*f*(loop*dt-100*dt))*1000*dz;
// Uz=-(sin(2*Pi*(loop-0.5)*f*dt)+sin(2*Pi*(loop-0.5)*(f+20000000)*dt)+sin(2*Pi*(loop-0.5)*(f+40000000)*dt)+sin(2*Pi*(loop-0.5)*(f+60000000)*dt)+sin(2*Pi*(loop-0.5)*(f+80000000)*dt)+sin(2*Pi*(loop-0.5)*(f+100000000)*dt)+sin(2*Pi*(loop-0.5)*(f+120000000)*dt)+sin(2*Pi*(loop-0.5)*(f+140000000)*dt)+sin(2*Pi*(loop-0.5)*(f+160000000)*dt)+sin(2*Pi*(loop-0.5)*(f+180000000)*dt)+sin(2*Pi*(loop-0.5)*(f+200000000)*dt)+sin(2*Pi*(loop-0.5)*(f+220000000)*dt)+sin(2*Pi*(loop-0.5)*(f+240000000)*dt)+sin(2*Pi*(loop-0.5)*(f+260000000)*dt)+sin(2*Pi*(loop-0.5)*(f+280000000)*dt)+sin(2*Pi*(loop-0.5)*(f+300000000)*dt))*dz*1000;
Uz1=-sin(2*Pi*f*(loop-0.5)*dt)*1000*dz;
Uzz1[loop-1]=Uz1;
Iz1=(dx*(Hxy[Nx/2][Ny/2-1][Nz/2]-Hxy[Nx/2][Ny/2][Nz/2])+dy*(Hyz[Nx/2][Ny/2][Nz/2]-Hyz[Nx/2-1][Ny/2][Nz/2]))*1000;
Izz1[loop-1]=Iz1;
Uz2=-sin(2*Pi*f*(loop-0.5)*dt)*1000*dz;
//Uz2=-sin(2*Pi*f*(loop-0.5)*dt)*1000*dz;
Uzz2[loop-1]=Uz2;
Iz2=(dy*(Hyz[Nx/2][Ny/2][Nz/2-1]-Hyz[Nx/2][Ny/2][Nz/2])+dz*(Hzx[Nx/2][Ny/2][Nz/2]-Hzx[Nx/2][Ny/2-1][Nz/2]))*1000;
Izz2[loop-1]=Iz2;
printf("%f,%f\n",Uz1,Iz1);
printf("%f,%f\n",Uz2,Iz2);
m=chuanbo2(loop);
printf("%d\n",loop);
printf("%f\n",Ezx[Nx/2][Ny/2][Nz/2]);
printf("%f\n",Ezy[Nx/2][Ny/2][Nz/2]);
// E1[loop]=sqrt(Ezx[Nx/2][Ny/2][Nz/2]*Ezx[Nx/2][Ny/2][Nz/2]+Exy[Nx/2][Ny/2][Nz/2]*Exy[Nx/2][Ny/2][Nz/2]+Eyz[Nx/2][Ny/2][Nz/2]*Eyz[Nx/2][Ny/2][Nz/2]);
// H1[loop]=sqrt(Hzx[Nx/2][Ny/2][Nz/2]*Hzx[Nx/2][Ny/2][Nz/2]+Hxy[Nx/2][Ny/2][Nz/2]*Hxy[Nx/2][Ny/2][Nz/2]+Hyz[Nx/2][Ny/2][Nz/2]*Hyz[Nx/2][Ny/2][Nz/2]);
flag3++;
if(loop>=NN-100)
{
for(loop1=0;loop1<=DD1+1;loop1++)
{
In1[loop1]=(dx*(Hxy[Nx/2][Ny/2-1][Nz/2+loop1-DD1/2-1]-Hxy[Nx/2][Ny/2][Nz/2+loop1-DD1/2-1])+dy*(Hyz[Nx/2][Ny/2][Nz/2+loop1-DD1/2-1]-Hyz[Nx/2-1][Ny/2][Nz/2+loop1-DD1/2-1]))*1000;
if(loop==NN-100)
{ current1[loop1].zhengfu=In1[loop1];current1[loop1].xiangwei=0;
}
else
{ if(In1[loop1]>current1[loop1].zhengfu)
{current1[loop1].zhengfu=In1[loop1];current1[loop1].xiangwei=(loop-NN+100)/100.0*2*Pi;}
}
}
for(loop1=0;loop1<=DD2+1;loop1++)
{
In2[loop1]=(dy*(Hyz[Nx/2+loop1-DD2/2-1][Ny/2+movedis][Nz/2-1]-Hyz[Nx/2+loop1-DD2/2-1][Ny/2+movedis][Nz/2])+dz*(Hzx[Nx/2+loop1-DD2/2-1][Ny/2+movedis][Nz/2]-Hzx[Nx/2+loop1-DD2/2-1][Ny/2-1+movedis][Nz/2]))*1000;
if(loop==NN-100)
{ current2[loop1].zhengfu=In2[loop1];current2[loop1].xiangwei=0;
}
else
{
if(In2[loop1]>current2[loop1].zhengfu)
{current2[loop1].zhengfu=In2[loop1];current2[loop1].xiangwei=(loop-NN+100)/100.0*2*Pi;}
}
}
}
if(loop>=NN-100)
distill(loop);
if(loop==NN)
waitui(loop);
}
printf("************************************\n");
printf("第%d 运算,第%d 次",DD1,DD2);
printf("************************************\n");
// ffft(Uzz1,Izz1);
// zukang();
zhubobi();
if ((out=fopen("jiexiao21","a"))==NULL)
{printf("cannot open outfile\n");
goto end;
}
else
{ fprintf(out,"两阵子长度为：%d z向，%d x向\n",DD1,DD2);
fprintf(out,"两阵子半径为：%f z向，%f x向\n",r0,r1);
for(j=0;j<=DD1+1;j++)
{
fprintf(out,"%f %f\n",current1[j].zhengfu,current1[j].xiangwei);
}
fprintf(out,"\n \n");
for(j=0;j<=DD2+1;j++)
fprintf(out,"%f %f\n",current2[j].zhengfu,current2[j].xiangwei);
}
if ((out=fopen("jiexiao22","a"))==NULL)
{printf("cannot open outfile\n");
goto end;
}
else
{ fprintf(out,"两阵子长度为：%d z向，%d x向\n",DD1,DD2);
fprintf(out,"两阵子半径为：%f z向，%f x向\n",(c/f)/r0,(c/f)/r1);
// fprintf(out,"阵子间距:%d",movedis);
fprintf(out,"%f %f;\n",fxhs.sita.shi,fxhs.sita.xu);
fprintf(out,"%f %f;\n",fxhs.fai.shi,fxhs.fai.xu);
fprintf(out,"\n\n\n");
}
}
end: return(0);
}
void fuzhi()
{int loop1,loop2,loop3;
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;
}
}
} /*给x向电场赋初值*/
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;
}
}
} /*给y向电场赋初值*/
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;
}
}
} /*给z向电场赋初值,并加激励源*/
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;
}
}
} /*给x向磁场赋初值*/
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;
}
}
} /*给y向磁场赋初值*/
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;
}
}
} /*给z向磁场赋初值*/
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;
}
fxhs.sita.shi=0;fxhs.sita.xu=0;
fxhs.fai.shi=0;fxhs.fai.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<Nx-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;
Jxl[loop1][loop2].zhengfu=0; Jxl[loop1][loop2].xiangwei=0;
Jxr[loop1][loop2].zhengfu=0; Jxr[loop1][loop2].xiangwei=0;
Jzl[loop1][loop2].zhengfu=0; Jzl[loop1][loop2].xiangwei=0;
Jzr[loop1][loop2].zhengfu=0; Jzr[loop1][loop2].xiangwei=0;
Jmxu[loop1][loop2].zhengfu=0; Jmxu[loop1][loop2].xiangwei=0;
Jmxd[loop1][loop2].zhengfu=0; Jmxd[loop1][loop2].xiangwei=0;
Jmxl[loop1][loop2].zhengfu=0; Jmxl[loop1][loop2].xiangwei=0;
Jmxr[loop1][loop2].zhengfu=0; Jmxr[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<Nx-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;
Jyf[loop1][loop2].zhengfu=0; Jyf[loop1][loop2].xiangwei=0;
Jyb[loop1][loop2].zhengfu=0; Jyb[loop1][loop2].xiangwei=0;
Jzf[loop1][loop2].zhengfu=0; Jzf[loop1][loop2].xiangwei=0;
Jzb[loop1][loop2].zhengfu=0; Jzb[loop1][loop2].xiangwei=0;
Jmyu[loop1][loop2].zhengfu=0; Jmyu[loop1][loop2].xiangwei=0;
Jmyd[loop1][loop2].zhengfu=0; Jmyd[loop1][loop2].xiangwei=0;
Jmyf[loop1][loop2].zhengfu=0; Jmyf[loop1][loop2].xiangwei=0;
Jmyb[loop1][loop2].zhengfu=0; Jmyb[loop1][loop2].xiangwei=0;
Jmzf[loop1][loop2].zhengfu=0; Jmzf[loop1][loop2].xiangwei=0;
Jmzb[loop1][loop2].zhengfu=0; Jmzb[loop1][loop2].xiangwei=0;
}
for(loop1=0;loop1<DD;loop1++)
{current1[loop1].zhengfu=0;current1[loop1].xiangwei=0;
current2[loop1].zhengfu=0;current2[loop1].xiangwei=0;}
for(loop1=0;loop1<DD;loop1++)
{In1[loop1]=0;In2[loop1]=0;}
}
/***************************************************************************************************************************************/
/*x向电场*/
int chuanbo1(int M) /*设置波动源，向前波动1个周期，并实现前时刻和后一时刻值的传递。M为时间点，需要返回的参数以便进行调试*/
{int loop1,loop2,loop3;
for(loop1=0;loop1<Nx;loop1++)
{for(loop2=1;loop2<Ny;loop2++)
{for(loop3=1;loop3<Nz;loop3++)
//if(loop1==N-NNN-5&&loop2==N-NNN-5&&loop3==N-NNN-5)
{ if((loop1==Nx/2||loop1==Nx/2-1)&&loop2==Ny/2+movedis&&loop3==Nz/2) //前一个5是预留网格,后一个5是介质厚度
{
Exy[loop1][loop2][loop3]=((2*e*er-sig*dt)/(2*e*er+sig*dt))*Exy[loop1][loop2][loop3]+(2*dt)/(2*e*er+sig*dt)*((Hzx[loop1][loop2][loop3]+Hzy[loop1][loop2][loop3]-Hzx[loop1][loop2-1][loop3]-Hzy[loop1][loop2-1][loop3])/dy)+sin(2*Pi*(M-0.5)*dt*f);
Exz[loop1][loop2][loop3]=((2*e*er-sig*dt)/(2*e*er+sig*dt))*Exz[loop1][loop2][loop3]-(2*dt)/(2*e*er+sig*dt)*((Hyx[loop1][loop2][loop3]+Hyz[loop1][loop2][loop3]-Hyz[loop1][loop2][loop3-1]-Hyx[loop1][loop2][loop3-1])/dz);
}
else
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -