antenna.c

贴片天线的FDTD三维编程

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define ie 100//62
#define je 100//120
#define ke 100//14
#define ia 7//14
#define ja 7//15
#define ka 7
#define ktop 3//2

main()
{
	    float dx[ie][je][ke],dy[ie][je][ke],dz[ie][je][ke];    
      	float ex[ie][je][ke],ey[ie][je][ke],ez[ie][je][ke];
	    float hx[ie][je][ke],hy[ie][je][ke],hz[ie][je][ke];
	    float gax[ie][je][ke],gay[ie][je][ke],gaz[ie][je][ke];

		int l=0,m=0,n=0,i=0,j=0,k=0,ic=0,jc=0,kc=0,nsteps=0,n_pml=0;
		float ddz=0,ra_x=0,ra_y=0,dt=0,T=0,epsz=0,muz=0,pi=0,eaf=0,npml=0;
		int ib=0,jb=0,kb=0;
		float xn=0,xxn=0,xnum=0,xd=0;
	    float t0=0,spread=0,pulse=0;
		float ez_low_m1=0,ez_low_m2=0,ez_high_m1=0,ez_high_m2=0;
		
        float ez_inc[je],hx_inc[je];
		float idxl[ia][je][ke],idxh[ia][je][ke];
		float ihxl[ia][je][ke],ihxh[ia][je][ke];
		float idyl[ie][ja][ke],idyh[ie][ja][ke];
		float ihyl[ie][ja][ke],ihyh[ie][ja][ke];
		float idzl[ie][je][ka],idzh[ie][je][ka];
		float ihzl[ie][je][ka],ihzh[ie][je][ka];
		int ixh=0,jyh=0,kzh=0;

		float gi1[ie],gi2[ie],gi3[ie];
        float gj1[je],gj2[je],gj3[je];
        float gk1[ke],gk2[ke],gk3[ke];
		float fi1[ie],fi2[ie],fi3[ie];
        float fj1[je],fj2[je],fj3[je];
        float fk1[ke],fk2[ke],fk3[ke];

		float curl_h=0,curl_e=0;

	    int istart=0,iend=0,k_ref=0,jend=0,i_ref=0;
        int j_patch_st=0,j_patch_end=0,j_ref=0;
        int i_patch_st=0,i_patch_end=0;
		float eps_sub=0;
		float shape[ie][ke],half_wv=0;

		FILE *fp,*fpt;

		ic=ie/2;
		jc=je/2;
		kc=ke/2;
		ib=ie-ia-1;
		jb=je-ja-1;
		kb=ke-ka-1;

		i_patch_st=ia+5;
		i_patch_end=i_patch_st+31;
		j_patch_end=jb-5;
		j_patch_st=j_patch_end-39;
		j_ref=j_patch_st-30;
		k_ref=ktop-1;

        pi=3.14159;
		epsz=8.8e-12;
		muz=4*pi*1.e-7;
		ddz=0.265e-3;             //cell size
		ra_y=0.6625;
		ra_x=0.6812;
		dt=ddz/6e8;               //time steps

		//initialize parameters calculate pml parameters
		//dielectric constant of the substrate material
		eps_sub=2.2;



		for(j=0;j<je;j++){
			for(i=0;i<ie;i++){
				for(k=0;k<ktop;k++){
				gax[i][j][k]=1.0/eps_sub;
				gay[i][j][k]=1.0/eps_sub;
				gaz[i][j][k]=1.0/eps_sub;
				}
			}
		}

		//add a metal plate at k=0
		for(i=1;i<ie-1;i++){
			for(j=1;j<je-1;j++){
			k=0;
			gax[i][j][k]=0.0;
			gay[i][j][k]=0.0;
			}
		}

		istart=ia+6;
		iend=istart+6;
		i_ref=istart+(iend-istart)/2;

		half_wv=(iend-istart)/2;

		for(i=istart;i<=iend;i++){
			for(k=0;k<=ktop;k++){
			shape[i][k]=1.0;
			}
		}

		//add the conductor lead at k=ktop+1
		for(j=1;j<=j_patch_st;j++){
			for(i=istart;i<=iend-1;i++){
			gax[i][j][ktop+1]=0.;
			gay[i][j][ktop+1]=0.;
			}
		}

		//add the rectangular patch at k=ktop
		for(j=j_patch_st;j<=j_patch_end;j++){
			for(i=ia+1;i<=ib-1;i++){
			gax[i][j][ktop+1]=0.;
			gay[i][j][ktop+1]=0.;
			}
		}

		t0=150.0;
		spread=25.0;
		T=0;
		nsteps=1;

		fpt=fopen("time.txt","w");

		while(nsteps>0){
			printf("nsteps-->");
			scanf("%d",&nsteps);
			
			for(n=1;n<=nsteps;n++){
			T=T+1;
			//start of the main fdtd loop
            //calculate the incident buffer
			for(j=1;j<je;j++){
			ez_inc[j]=gj3[j]*ez_inc[j]+gj2[j]*(0.5*ra_y/eps_sub)*(hx_inc[j-1]-hx_inc[j]);
			}

           //source
			pulse=exp(-0.5*(pow((t0-T)/spread,2.0)));
			ez_inc[ja-4]=pulse;

			//calculate the Dx field
            for(i=1;i<ia;i++){            
				for(j=1;j<je;j++){
					for(k=1;k<ke;k++){
					curl_h=(ra_y*(hz[i][j][k]-hz[i][j-1][k])-hy[i][j][k]+hy[i][j][k-1]);
					idxl[i][j][k]=idxl[i][j][k]+curl_h;
                    dx[i][j][k]=gj3[j]*gk3[k]*dx[i][j][k]+gj2[j]*gk2[k]*0.5*(curl_h+gi1[i]*idxl[i][j][k]);
					}
				}
			}

			for(i=ia;i<ib;i++){            
				for(j=1;j<je;j++){
					for(k=1;k<ke;k++){
					curl_h=(ra_y*(hz[i][j][k]-hz[i][j-1][k])-hy[i][j][k]+hy[i][j][k-1]);
                    dx[i][j][k]=gj3[j]*gk3[k]*dx[i][j][k]+gj2[j]*gk2[k]*0.5*curl_h;
					}
				}
			}

		 for(i=ib+1;i<ie;i++){   
			 ixh=i-ib-1;
				for(j=1;j<je;j++){
					for(k=1;k<ke;k++){
					curl_h=(ra_y*(hz[i][j][k]-hz[i][j-1][k])-hy[i][j][k]+hy[i][j][k-1]);
					idxh[i][j][k]=idxh[i][j][k]+curl_h;
                    dx[i][j][k]=gj3[j]*gk3[k]*dx[i][j][k]+gj2[j]*gk2[k]*0.5*(curl_h+gi1[i]*idxh[ixh][j][k]);
					}
				}
			}

		 //claculate the Dy field
            for(i=1;i<ie;i++){            
				for(j=1;j<ja;j++){
					for(k=1;k<ke;k++){
					curl_h=(hx[i][j][k]-hx[i][j][k-1]-ra_x*(hz[i][j][k]-hz[i-1][j][k]));
					idyl[i][j][k]=idyl[i][j][k]+curl_h;
                    dy[i][j][k]=gi3[i]*gk3[k]*dy[i][j][k]+gi2[i]*gk2[k]*0.5*(curl_h+gj1[j]*idyl[i][j][k]);
					}
				}
			}

            for(i=1;i<ie;i++){            
				for(j=ja;j<=jb;j++){
					for(k=1;k<ke;k++){
					curl_h=(hx[i][j][k]-hx[i][j][k-1]-ra_x*(hz[i][j][k]-hz[i-1][j][k]));
				    dy[i][j][k]=gi3[i]*gk3[k]*dy[i][j][k]+gi2[i]*gk2[k]*0.5*curl_h;
					}
				}
			}
			
		for(i=1;i<ie;i++){            
				for(j=jb+1;j<je;j++){
					jyh=j-jb-1;
					for(k=1;k<ke;k++){
					curl_h=(hx[i][j][k]-hx[i][j][k-1]-ra_x*(hz[i][j][k]-hz[i-1][j][k]));
					idyh[i][jyh][k]=idyh[i][jyh][k]+curl_h;
                    dy[i][j][k]=gi3[i]*gk3[k]*dy[i][j][k]+gi2[i]*gk2[k]*0.5*(curl_h+gj1[j]*idyh[i][jyh][k]);
					}
				}
			}

//calculate the Dz field
		     for(i=1;i<ie;i++){            
				for(j=1;j<je;j++){
					for(k=1;k<ka;k++){
					curl_h=(ra_x*(hy[i][j][k]-hy[i-1][j][k])-ra_y*(hx[i][j][k]-hx[i][j-1][k]));
					idzl[i][j][k]=idzl[i][j][k]+curl_h;
                    dz[i][j][k]=gi3[i]*gj3[j]*dz[i][j][k]+gi2[i]*gj2[j]*0.5*(curl_h+gk1[k]*idzl[i][j][k]);
					}
				}
			}

		     for(i=1;i<ie;i++){            
				for(j=1;j<je;j++){
					for(k=ka;k<kb;k++){
					curl_h=(ra_x*(hy[i][j][k]-hy[i-1][j][k])-ra_y*(hx[i][j][k]-hx[i][j-1][k]));
				    dz[i][j][k]=gi3[i]*gj3[j]*dz[i][j][k]+gi2[i]*gj2[j]*0.5*curl_h;
					}
				}
			}

		     for(i=1;i<ie;i++){            
				for(j=1;j<je;j++){
					for(k=kb+1;k<ke;k++){
					kzh=k-kb-1;
					curl_h=(ra_x*(hy[i][j][k]-hy[i-1][j][k])-ra_y*(hx[i][j][k]-hx[i][j-1][k]));
					idzh[i][j][kzh]=idzh[i][j][kzh]+curl_h;
                    dz[i][j][k]=gi3[i]*gj3[j]*dz[i][j][k]+gi2[i]*gj2[j]*0.5*(curl_h+gk1[k]*idzh[i][j][kzh]);
					}
				}
			}

			 //incident Dz
			 for(i=istart;i<=iend;i++){
				 for(k=0;k<=ktop;k++){
				 dz[i][ja][k]=dz[i][ja][k]+(0.5/eps_sub)*shape[i][k]*hx_inc[ja-1];
				 }
			 }

			 //calculate the E form D field
			 //Ex and Ey are zero at k=0.
			 for(i=1;i<ie-1;i++){
				 for(j=1;j<je-1;j++){
					 for(k=1;k<ke-1;k++){
					 ex[i][j][k]=gax[i][j][k]*dx[i][j][k];
					 ey[i][j][k]=gay[i][j][k]*dy[i][j][k];
					 ez[i][j][k]=gaz[i][j][k]*dz[i][j][k];
					 }
				 }
			 }

			 //write out the time domain data at the input port
			 //clculate the incident field
			 for(j=0;j<je-1;j++){
			 hx_inc[j]=fj3[j]*hx_inc[j]+0.5*fj2[j]*(ez_inc[j]-ez_inc[j+1]);
			 }

			 //calculate the Hx field
			 for(i=0;i<ia;i++){
				 for(j=0;j<je-1;j++){
					 for(k=0;k<ke-1;k++){
				 curl_e=(ey[i][j][k+1]-ey[i][j][k]-ra_y*(ez[i][j+1][k]-ez[i][j][k]));
				 ihxl[i][j][k]=ihxl[i][j][k]+curl_e;
				 hx[i][j][k]=fj3[j]*fk3[k]*hx[i][j][k]+fj2[j]*fk2[k]*0.5*(curl_e+fi1[i]*ihxl[i][j][k]);
					 }
				 }
			 }

			 for(i=ia;i<=ib;i++){
				 for(j=0;j<je-1;j++){
					 for(k=0;k<ke-1;k++){
				 curl_e=(ey[i][j][k+1]-ey[i][j][k]-ra_y*(ez[i][j+1][k]-ez[i][j][k]));
				 ihxl[i][j][k]=ihxl[i][j][k]+curl_e;
				 hx[i][j][k]=fj3[j]*fk3[k]*hx[i][j][k]+fj2[j]*fk2[k]*0.5*curl_e;
					 }
				 }
			 }
			 for(i=ia;i<=ib;i++){
				 for(j=0;j<je-1;j++){
					 for(k=0;k<ke-1;k++){
				 curl_e=(ey[i][j][k+1]-ey[i][j][k]-ra_y*(ez[i][j+1][k]-ez[i][j][k]));
				 hx[i][j][k]=fj3[j]*fk3[k]*hx[i][j][k]+fj2[j]*fk2[k]*0.5*curl_e;
					 }
				 }
			 }

			 for(i=ib+1;i<ie;i++){
				 ixh=i-ib-1;
				 for(j=0;j<je-1;j++){
					 for(k=0;k<ke-1;k++){
				 curl_e=(ey[i][j][k+1]-ey[i][j][k]-ra_y*(ez[i][j+1][k]-ez[i][j][k]));
				 ihxh[ixh][j][k]=ihxh[ixh][j][k]+curl_e;
				 hx[i][j][k]=fj3[j]*fk3[k]*hx[i][j][k]+fj2[j]*fk2[k]*0.5*(curl_e+fi1[i]*ihxh[ixh][j][k]);
					 }
				 }
			 }

			 //incident Hx
			 for(i=istart;i<=iend;i++){
				 for(k=0;k<=ktop;k++){
				 hx[i][ja-1][k]=hx[i][ja-1][k]+(0.5/eps_sub)*shape[i][k]*ez_inc[ja];
				 }
			 }


			 //calculate the Hy field
			 for(i=0;i<ie-1;i++){
				 for(j=0;j<ja;j++){
					 for(k=0;k<ke-1;k++){
					 curl_e=(ra_x*(ez[i+1][j][k]-ez[i][j][k])-ex[i][j][k+1]+ex[i][j][k]);
					 ihyl[i][j][k]=ihyl[i][j][k]+curl_e;
                     hy[i][j][k]=fi3[i]*fk3[k]*hy[i][j][k]+fi2[i]*fk2[k]*0.5*(curl_e+fj1[j]*ihyl[i][j][k]);
					 }
				 }
			 }

			 for(i=0;i<ie-1;i++){
				 for(j=ja;j<=jb;j++){
					 for(k=0;k<ke-1;k++){
					 curl_e=(ra_x*(ez[i+1][j][k]-ez[i][j][k])-ex[i][j][k+1]+ex[i][j][k]);
                     hy[i][j][k]=fi3[i]*fk3[k]*hy[i][j][k]+fi2[i]*fk3[k]*0.5*curl_e;
					 }
				 }
			 }


			 for(i=0;i<ie-1;i++){
				 for(j=jb+1;j<je;j++){
					 jyh=j-jb-1;
					 for(k=0;k<ke-1;k++){
					 curl_e=(ra_x*(ez[i+1][j][k]-ez[i][j][k])-ex[i][j][k+1]+ex[i][j][k]);
					 ihyh[i][jyh][k]=ihyh[i][jyh][k]+curl_e;
                     hy[i][j][k]=fi3[i]*fk3[k]*hy[i][j][k]+fi2[i]*fk2[k]*0.5*(curl_e+fj1[j]*ihyh[i][jyh][k]);
					 }
				 }
			 }

			 //calculate the Hz field
			 for(i=0;i<ie-1;i++){
				 for(j=0;j<je-1;j++){
					 for(k=0;k<ka;k++){
					 curl_e=(ra_y*(ex[i][j+1][k]-ex[i][j][k])-ra_x*(ey[i+1][j][k]-ey[i][j][k]));
					 ihzl[i][j][k]=ihzl[i][j][k]+curl_e;
                     hz[i][j][k]=fi3[i]*fj3[j]*hz[i][j][k]+fi2[i]*fj2[j]*0.5*(curl_e+fk1[k]*ihzl[i][j][k]);
					 }
				 }
			 }

			 for(i=0;i<ie-1;i++){
				 for(j=0;j<je-1;j++){
					 for(k=ka;k<=kb;k++){
					 curl_e=(ra_y*(ex[i][j+1][k]-ex[i][j][k])-ra_x*(ey[i+1][j][k]-ey[i][j][k]));
                     hz[i][j][k]=fi3[i]*fj3[j]*hz[i][j][k]+fi2[i]*fj2[j]*0.5*curl_e;
					 }
				 }
			 }

			 for(i=0;i<ie-1;i++){
				 for(j=0;j<je-1;j++){
					 for(k=kb+1;k<ke;k++){
					kzh=k-kb-1;
					 curl_e=(ra_y*(ex[i][j+1][k]-ex[i][j][k])-ra_x*(ey[i+1][j][k]-ey[i][j][k]));
					 ihzl[i][j][k]=ihzl[i][j][k]+curl_e;
                     hz[i][j][k]=fi3[i]*fj3[j]*hz[i][j][k]+fi2[i]*fj2[j]*0.5*(curl_e+fk1[k]*ihzl[i][j][k]);
					 }
				 }
			 }
}
			 //end of the main fdtd loop
			 //write the E field out to a file "Ez"
             fp=fopen("Ez.txt","w");
			 for(j=0;j<je;j++){
				 for(i=0;i<ie;i++){
				 fprintf(fp,"%9.6f",ez[i][j][ktop]);
				 }
				 fprintf(fp,"\n");
			 }

              fclose(fp);

			}

         fclose(fpt);

		}