?? xiangyang28.cpp
字號:
/*PML完全匹配層*/
#include "shuju9.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));
DD2=50;
for(DD1=16;DD1<24;DD1=DD1+2)
{
fuzhi();
// for (movedis=-10;movedis<=10;movedis=movedis+5)
// {
// fuzhi();
for(lp=0;lp<10;lp++)
{ r0=(c/f)/500.0;r1=(c/f)/(100*(lp+1));
// 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));
for(loop=1;loop<=NN;loop++) /*主時間循環(huán)*/
{
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);
if(r1==(c/f)/(100)||r1==(c/f)/(200)||r1==(c/f)/(300))
{m=chuanbo3(loop);}
else
{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("jiexiao3","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("jiexiao4","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<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;
Jmyf[loop1][loop2].zhengfu=0; Jmyf[loop1][loop2].xiangwei=0;
Jmyb[loop1][loop2].zhengfu=0; Jmyb[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++)
{
Jzf[loop1][loop2].zhengfu=0; Jzf[loop1][loop2].xiangwei=0;
Jzb[loop1][loop2].zhengfu=0; Jzb[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個周期，并實現(xiàn)前時刻和后一時刻值的傳遞。M為時間點，需要返回的參數(shù)以便進行調試*/
{int loop1,loop2,loop3;
for(loop1=0;loop1<Nx;loop1++)
{for(loop2=1;loop2<Ny;loop2++)
{for(loop3=1;loop3<Nz;loop3++)
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -