📄 patch1.cpp
字号:
for (k=0;k<(KE-1);k++){
Hxy[i][j][k]=CPy[jy]*Hxy[i][j][k]-CQy[jy]*(Ez[i][j+1][k]-Ez[i][j][k]);
Hx[i][j][k]=Hxy[i][j][k]+Hxz[i][j][k];
}}
for (k=1;k<(KE-1);k++){
for (i=n_pml;i<i_max;i++){
Hzx[i][j][k]=Hzx[i][j][k]-.5*ra_x*(Ey[i+1][j][k]-Ey[i][j][k]);
}
for (i=0;i<(IE-1);i++){
Hzy[i][j][k]=CPy[jy]*Hzy[i][j][k]+CQy[jy]*(Ex[i][j+1][k]-Ex[i][j][k]);
Hz[i][j][k]=Hzx[i][j][k]+Hzy[i][j][k];
}}}
/* j=j_max:(JE-1) */
for (j=j_max;j<(JE-1);j++){
jy=j-j_max;
for (i=1;i<(IE-1);i++){
for (k=0;k<k_max;k++){
Hxz[i][j][k]=Hxz[i][j][k]+.5*(Ey[i][j][k+1]-Ey[i][j][k]);
}
for (k=0;k<(KE-1);k++){
Hxy[i][j][k]=CPy[jy]*Hxy[i][j][k]-CQy[jy]*(Ez[i][j+1][k]-Ez[i][j][k]);
Hx[i][j][k]=Hxy[i][j][k]+Hxz[i][j][k];
}}
for (k=1;k<(KE-1);k++){
for (i=n_pml;i<i_max;i++){
Hzx[i][j][k]=Hzx[i][j][k]-.5*ra_x*(Ey[i+1][j][k]-Ey[i][j][k]);
}
for (i=0;i<(IE-1);i++){
Hzy[i][j][k]=CPy[jy]*Hzy[i][j][k]+CQy[jy]*(Ex[i][j+1][k]-Ex[i][j][k]);
Hz[i][j][k]=Hzx[i][j][k]+Hzy[i][j][k];
}}}
/* Write out the time domain data at J_Jm input port */
/*printf("%4.0d %6.2f \n",n,Ez[i_ref][j_ref][k_ref]);*/
printf("%4d ",n);
fprintf(fp,"%4.0d %f\n",n,Ez[i_ref][j_ref][k_ref]);
/* Store the imaginary part of V,I,J and Jm */
if (n==(nsteps-int(1./f0/dt/4.+.5))){
aaa=0;
for (k=0;k<k_top;k++){
aaa=aaa+ddz*Ez[i_ref][j_ref][k];
}
Voltage=Complex(0,aaa);
aaa=0;
for (i=i_mic_start;i<i_mic_end;i++){
aaa=aaa+ddx*Hx[i][j_ref][k_top-1]-ddx*Hx[i][j_ref][k_top];
}
aaa=aaa-ddz*Hz[i_mic_start-1][j_ref][k_top]+ddz*Hz[i_mic_end-1][j_ref][k_top];
Current=Complex(0,aaa);
for (i=0;i<IE;i++){
for (j=0;j<JE;j++){
Jx[i][j]=-double(Hy[i][j][k_top])*_J;
Jy[i][j]=double(Hx[i][j][k_top])*_J;
Jmx[i][j]=double(Ey[i][j][k_top])*_J;
Jmy[i][j]=-double(Ex[i][j][k_top])*_J;
}}
}
}
fclose (fp);
/* End of the main FDTD loop */
/*Write the Efield out to a tile "Ez" */
fp=fopen(Ez_file,"w");
for (j=0;j<JE;j++){
for (i=0;i<IE;i++){
fprintf(fp,"%f ",Ez[i][j][k_top-1]);
}
fprintf(fp," \n");
}
fclose(fp);
/* Store the real part of V,I,J and Jm */
aaa=0;
for (k=0;k<k_top;k++){
aaa=aaa+ddz*Ez[i_ref][j_ref][k];
}
Voltage=Voltage+Complex(aaa,0);
aaa=0;
for (i=i_mic_start;i<i_mic_end;i++){
aaa=aaa+ddx*Hx[i][j_ref][k_top-1]-ddx*Hx[i][j_ref][k_top];
}
aaa=aaa-ddz*Hz[i_mic_start-1][j_ref][k_top]+ddz*Hz[i_mic_end-1][j_ref][k_top];
Current=Current+Complex(aaa,0);
for (j=0;j<JE;j++){
for (i=0;i<IE;i++){
Jx[i][j]=Jx[i][j]-double(Hy[i][j][k_top]);
Jy[i][j]=Jy[i][j]+double(Hx[i][j][k_top]);
Jmx[i][j]=Jmx[i][j]+double(Ey[i][j][k_top]);
Jmy[i][j]=Jmy[i][j]-double(Ex[i][j][k_top]);
}
}
Zin=Complex(376.8194)*Voltage/Current;
S11=(Z0-Zin)/(Z0+Zin);
cal_pattern();
cout<<'\n';
printf("IE:%d\; JE:%d\; KE:%d\; \n",IE,JE,KE);
printf(Ez_file);
printf("\n");
printf(time_file);
printf("\n");
cout<<"Zin:"<<Zin<<" S11:"<<S11<<" "<<abs(S11)<<" "<<20*log10(abs(S11))<<"(dB)"<<'\n';
outfile.open(S11_file,ios::out);
outfile<<"Zin:"<<Zin<<" S11:"<<S11<<" "<<abs(S11)<<" "<<20*log10(abs(S11))<<"(dB)"<<'\n';
outfile.close();
time_finish = clock();
duration = (double)(time_finish - time_start)/CLOCKS_PER_SEC;
printf( "\nTime to finish this FDTD program is %f seconds\n", duration);
}
void cal_pattern()
{ /* calculate the pattern */
double ds;
double E_avr;
int t,p;
double theta,phi;
ds=(ddx*ddy);
E_avr=sqrt(real(Complex(real(Current),-imag(Current))*Voltage)*376.8194/4/pi*376.8194);
outfile.open(Pattern_file,ios::out);
for (t=0;t<=180;t++){
theta=pi/2.-t*pi/180;
for (p=0;p<=1;p++){
phi=p*90*pi/180.;
fx[t][p]=Complex(0,0);
fy[t][p]=Complex(0,0);
fmx[t][p]=Complex(0,0);
fmy[t][p]=Complex(0,0);
for (i=(n_pml);i<=i_max;i++){
for (j=j_ref;j<=j_max;j++){
fx[t][p]=fx[t][p]+Jx[i][j]*ds*exp(_J*(k0*(i-i_center+.5)*ddx*sin(theta)*cos(phi)+k0*(j-j_center)*ddy*sin(theta)*sin(phi)));
fy[t][p]=fy[t][p]+Jy[i][j]*ds*exp(_J*(k0*(i-i_center)*ddx*sin(theta)*cos(phi)+k0*(j-j_center+.5)*ddy*sin(theta)*sin(phi)));
fmx[t][p]=fmx[t][p]+Jmx[i][j]*ds*exp(_J*(k0*(i-i_center)*ddx*sin(theta)*cos(phi)+k0*(j-j_center+.5)*ddy*sin(theta)*sin(phi)));
fmy[t][p]=fmy[t][p]+Jmy[i][j]*ds*exp(_J*(k0*(i-i_center+0.5)*ddx*sin(theta)*cos(phi)+k0*(j-j_center)*ddy*sin(theta)*sin(phi)));
}}
fx[t][p]=fx[t][p]*(376.8194*2.*sin((k_top+.5)*ddz*k0*cos(theta))*_J);
fy[t][p]=fy[t][p]*(376.8194*2.*sin((k_top+.5)*ddz*k0*cos(theta))*_J);
fmx[t][p]=fmx[t][p]*(2.*cos(k_top*ddz*k0*cos(theta)));
fmy[t][p]=fmy[t][p]*(2.*cos(k_top*ddz*k0*cos(theta)));
fmx[t][p]=376.8194*fmx[t][p]*exp(_J*(2*pi*f0*dt/2));
fmy[t][p]=376.8194*fmy[t][p]*exp(_J*(2*pi*f0*dt/2));
E_theta[t][p]=-_J*k0*exp(-_J*k0)/(4*pi)*((fx[t][p]*(cos(theta)*cos(phi))+fy[t][p]*(cos(theta)*sin(phi)))+(fmx[t][p]*(-sin(phi))+fmy[t][p]*(cos(phi))));
E_phi[t][p]=_J*k0*exp(-_J*k0)/(4*pi)*((fx[t][p]*sin(phi)-fy[t][p]*cos(phi))+(fmx[t][p]*(cos(theta)*cos(phi))+fmy[t][p]*(cos(theta)*sin(phi))));
E_theta[t][p]=E_theta[t][p]/E_avr;
E_phi[t][p]=E_phi[t][p]/E_avr;
E_total[t][p]=sqrt(abs(E_theta[t][p]*E_theta[t][p])+abs(E_phi[t][p]*E_phi[t][p]));
outfile<<20*log10(abs(E_theta[t][p]))<<" "<<20*log10(abs(E_phi[t][p]))<<" "<<20*log10(abs(E_total[t][p]))<<" ";
}
outfile<<'\n';
}
outfile.close();
}
void init_parameter(){
if ((ra_x>1)||(ra_y>1)){
printf("dx or dy can not be larger than dz!\n");
exit(1);
}
strcat(Ez_file,This_file);
strcat(Pattern_file,This_file);
strcat(time_file,This_file);
strcat(S11_file,This_file);
/* init the parameter */
/*Dielectric constant of the substrate materi_minl */
for(j=0;j<JE;j++){
for(i=0;i<IE;i++){
k=0;
gaz[i][j][k]=1./sub_eps;
gax[i][j][k]=0;
gay[i][j][k]=0;
for(k=1;k<=k_top-1;k++){
gax[i][j][k]=1./sub_eps;
gay[i][j][k]=1./sub_eps;
gaz[i][j][k]=1./sub_eps;
}
k=k_top;
gax[i][j][k]=2./(sub_eps+1);
gay[i][j][k]=2./(sub_eps+1);
gaz[i][j][k]=1.;
for(k=k_top+1;k<KE;k++){
gax[i][j][k]=1.;
gay[i][j][k]=1.;
gaz[i][j][k]=1.;
}
}}
/* Initialize the arrays */
for (i=0;i<IE;i++){
for (j=0;j<JE;j++){
for (k=0;k<KE;k++){
Ex[i][j][k]=0.0;
Ey[i][j][k]=0.0;
Ez[i][j][k]=0.0;
Dx[i][j][k]=0.0;
Dy[i][j][k]=0.0;
Dz[i][j][k]=0.0;
Hx[i][j][k]=0.0;
Hy[i][j][k]=0.0;
Hz[i][j][k]=0.0;
Dxy[i][j][k]=0.0;
Dxz[i][j][k]=0.0;
Dyz[i][j][k]=0.0;
Dyx[i][j][k]=0.0;
Dzx[i][j][k]=0.0;
Dzy[i][j][k]=0.0;
Hxy[i][j][k]=0.0;
Hxz[i][j][k]=0.0;
Hyz[i][j][k]=0.0;
Hyx[i][j][k]=0.0;
Hzx[i][j][k]=0.0;
Hzy[i][j][k]=0.0;
}}}
/* Boundary conditions (PML)*/
float n=2;
float npml=n_pml,eps;
float L;
float oz,omz,sigma_z0,ox,omx,sigma_x0,oy,omy,sigma_y0;
/* z dirrection */
eps=eps0;
sigma_z0=eps*cc*(-log(R0))/(pow(2.,(n+2.))*ddz*pow(npml,(n+1)));
k=0;
oz=sigma_z0;
CAz[k]=exp(-oz*dt/eps);
CBz[k]=0.5*exp(-oz*dt/(2*eps));
for (k=1;k<n_pml;k++){
L=k;
oz=sigma_z0*(pow((2*L+1),(n+1))-pow((2*L-1),(n+1)));
CAz[k]=exp(-oz*dt/eps);
CBz[k]=0.5*exp(-oz*dt/(2*eps));
}
for (k=0;k<n_pml;k++){
L=k+0.5;
oz=sigma_z0*(pow((2*L+1),(n+1))-pow((2*L-1),(n+1)));
omz=oz/eps*mu0;
CPz[k]=exp(-omz*dt/mu0);
CQz[k]=.5*exp(-omz*dt/(2*mu0));
}
/* x direction */
eps=eps0;
sigma_x0=eps*cc*(-log(R0))/(pow(2.,(n+2.))*ddx*pow(npml,(n+1.)));
i=0;
ox=sigma_x0;
CAx[i]=exp(-ox*dt/eps);
CBx[i]=0.5*ra_x*exp(-ox*dt/(2*eps));
for (i=1;i<n_pml;i++){
L=i;
ox=sigma_x0*(pow((2*L+1),(n+1))-pow((2*L-1),(n+1)));
CAx[i]=exp(-ox*dt/eps);
CBx[i]=0.5*ra_x*exp(-ox*dt/(2*eps));
}
for (i=0;i<n_pml;i++){
L=i+0.5;
ox=sigma_x0*(pow((2*L+1),(n+1))-pow((2*L-1),(n+1)));
omx=ox/eps*mu0;
CPx[i]=exp(-omx*dt/mu0);
CQx[i]=.5*ra_x*exp(-omx*dt/(2*mu0));
}
/* y direction */
eps=eps0;
sigma_y0=eps*cc*(-log(R0))/(pow(2.,(n+2))*ddy*pow(npml,(n+1)));
j=0;
oy=sigma_y0;
CAy[j]=exp(-oy*dt/eps);
CBy[j]=0.5*ra_y*exp(-oy*dt/(2*eps));
for (j=1;j<n_pml;j++){
L=j;
oy=sigma_y0*(pow((2*L+1),(n+1))-pow((2*L-1),(n+1)));
CAy[j]=exp(-oy*dt/eps);
CBy[j]=0.5*ra_y*exp(-oy*dt/(2*eps));
}
for (j=0;j<n_pml;j++){
L=j+0.5;
oy=sigma_y0*(pow((2*L+1),(n+1))-pow((2*L-1),(n+1)));
omy=oy/eps*mu0;
CPy[j]=exp(-omy*dt/mu0);
CQy[j]=.5*ra_y*exp(-omy*dt/(2*mu0));
}
}
void show_parameter()
{
for (k=0;k<n_pml;k++){
cout<<"CAz:"<<CAz[k]<<" CBz:"<<CBz[k];
cout<<" CPz:"<<CPz[k]<<" CQz:"<<CQz[k]<<'\n';
}
system("pause");
for (k=0;k<KE;k++){
cout<<"k="<<k<<'\n';
for (i=0;i<n_pml;i++){
cout<<"CAx:"<<CAx[i]<<" CBx:"<<CBx[i];
cout<<" CPx:"<<CPx[i]<<" CQx:"<<CQx[i]<<'\n';
}
system("pause");
}
for (k=0;k<KE;k++){
cout<<"k="<<k<<'\n';
for (j=0;j<n_pml;j++){
cout<<"CAy:"<<CAy[j]<<" CBy:"<<CBy[j];
cout<<" CPy:"<<CPy[j]<<" CQy:"<<CQy[j]<<'\n';
}
system("pause");
}
}
void init_conductor(){
/* Add the conductor lead at k=k_top */
/* the microstrip conductor plane */
for (j=0;j<j_patch_st;j++){
for (i=i_mic_start;i<i_mic_end;i++){
gax[i][j][k_top]=0;
}}
for (j=0;j<j_patch_st;j++){
for (i=i_mic_start;i<i_mic_end;i++){
gay[i][j][k_top]=0;
}}
for (j=j_patch_st;j<j_patch_end;j++){
for (i=i_patch_st;i<i_patch_end;i++){
gax[i][j][k_top]=0;
}}
for (j=j_patch_st;j<j_patch_end;j++){
for (i=i_patch_st;i<i_patch_end;i++){
gay[i][j][k_top]=0;
}}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -