patch_s11.cpp

一个用C＋＋计算微带天线的完整的3D例子

		/* i=i_max:(IE-1) */
		for (i=i_max;i<(IE-1);i++){
			ix=i-i_max;
			for (j=1;j<(JE-1);j++){
				for (k=0;k<k_max;k++){
					Hyz[i][j][k]=Hyz[i][j][k]-.5*(Ex[i][j][k+1]-Ex[i][j][k]);
				}
				for (k=0;k<(KE-1);k++){
					Hyx[i][j][k]=CPx[ix]*Hyx[i][j][k]+CQx[ix]*(Ez[i+1][j][k]-Ez[i][j][k]);
					Hy[i][j][k]=Hyz[i][j][k]+Hyx[i][j][k];
			}}
			for (k=1;k<(KE-1);k++){
				for (j=0;j<(JE-1);j++){
					Hzx[i][j][k]=CPx[ix]*Hzx[i][j][k]-CQx[ix]*(Ey[i+1][j][k]-Ey[i][j][k]);
				}
				for (j=n_pml;j<j_max;j++){
					Hzy[i][j][k]=Hzy[i][j][k]+.5*ra_y*(Ex[i][j+1][k]-Ex[i][j][k]);
					Hz[i][j][k]=Hzx[i][j][k]+Hzy[i][j][k];
		}}}
		/* j=0:j_min-1 */
		for (j=0;j<n_pml;j++){
			jy=n_pml-j-1;
			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];
		}}}

		/* 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 eh 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]);

		V_all[n-1]=0;
		I_all[n-1]=0;
		for (k=0;k<k_top;k++)
		{	V_all[n-1]=V_all[n-1]+ddz*Ez[i_ref][j_ref][k];	}
		for (i=i_mic_start;i<i_mic_end;i++){
			I_all[n-1]=I_all[n-1]+ddx*Hx[i][j_ref][k_top-1]-ddx*Hx[i][j_ref][k_top];
		}
		I_all[n-1]=I_all[n-1]-ddz*Hz[i_mic_start-1][j_ref][k_top]+ddz*Hz[i_mic_end-1][j_ref][k_top];

	}
	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 */

	cal_s11();
	
	cout<<'\n';
	printf("IE:%d\; JE:%d\; KE:%d\; \n",IE,JE,KE);

	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_s11()
{	/* calculate the S11 */
	Complex V_in[f_count],V_out[f_count];
	Complex Voltage[f_count],Current[f_count];
	Complex S11[f_count],Zin[f_count];
	float freq[f_count],abs_S11[f_count];
	float f_step;

	f_step=(f_max-f_min)/(f_count-1);
	for (i=0;i<f_count;i++){
		V_in[i]=Complex(0.,0.);
		V_out[i]=Complex(0.,0.);
		Voltage[i]=Complex(0.,0.);
		Current[i]=Complex(0.,0.);
		freq[i]=f_min+i*f_step;
	}

	for (i=0;i<f_count;i++){
		for (n=0;n<n_cut;n++){
			V_in[i]=V_in[i]+V_all[n]*exp(Complex(0,-(2*pi*freq[i]*dt*n)));
			Current[i]=Current[i]+I_all[n]*exp(Complex(0,-(2*pi*freq[i]*dt*n)));
		}
		for (n=n_cut;n<nsteps;n++){
			V_out[i]=V_out[i]+V_all[n]*exp(Complex(0,-(2*pi*freq[i]*dt*n)));
			Current[i]=Current[i]+I_all[n]*exp(Complex(0,-(2*pi*freq[i]*dt*n)));
		}
		Voltage[i]=V_in[i]+V_out[i];
	}
	fp=fopen(S11_file,"w");
	for (i=0;i<f_count;i++){
		S11[i]=V_out[i]/V_in[i];
		abs_S11[i]=abs(S11[i]);
		Zin[i]=Voltage[i]/Current[i]*Complex(376.8194,0);
		fprintf(fp,"%f  %f  %f  %f  %f  %f\n",freq[i],real(Zin[i]),imag(Zin[i]),real(S11[i]),imag(S11[i]),abs_S11[i]);
	}
	fclose(fp);
}

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(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;
	}}
}