📄 main.c
字号:
/* fd3d_4.3c 3D Radiation from a dipole an FDTD program with
a senven_point PML and total/scantteried field formulation*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define IE 100 /*IE is the number of cells to the X lable*/
#define JE 140 /*IE is the number of cells to the Y lable*/
#define KE 100 /*IE is the number of cells to the Z lable*/
#define ia 15
#define ja 15
#define ka 15
#define pi 3.1415926
#define epsz 8.85419e-12
#define TT 500
#define LL 7
main()
{
static float dx[IE][JE][KE],dy[IE][JE][KE],dz[IE][JE][KE];
static float ex[IE][JE][KE],ey[IE][JE][KE],ez[IE][JE][KE];
static float hx[IE][JE][KE],hy[IE][JE][KE],hz[IE][JE][KE];
static float gax[IE][JE][KE],gay[IE][JE][KE],gaz[IE][JE][KE];
static int i,j,n,jc,ic,kc,k,NSTEPS,npml,n_pml,jmetal;
static int ib,jb,kb,kzh,ixh,jyh,T;
static float xn,xxn,curl_h,curl_e,curl_d;
static float ddx,dt,epsilon,omega;
static float t0,spread,pulse;
FILE *fp,*fopen(),*fptime,*fpt1,*fpt2,*fpt3,*fpt4,*fpt5,*fpt6,*fpt7;
static float ez_inc[JE],hx_inc[JE];
static float ez_low_m1=0.0,ez_low_m2=0.0,ez_high_m1=0.0,ez_high_m2=0.0;
static float gi1[IE],gi2[IE],gi3[IE];
static float gj1[JE],gj2[JE],gj3[JE];
static float gk1[KE],gk2[KE],gk3[KE];
static float idxl[ia][JE][KE],idxh[ia][JE][KE];
static float ihxl[ia][JE][KE],ihxh[ia][JE][KE];
static float idyl[IE][ja][KE],idyh[IE][ja][KE];
static float ihyl[IE][ja][KE],ihyh[IE][ja][KE];
static float idzl[IE][JE][ka],idzh[IE][JE][ka];
static float ihzl[IE][JE][ka],ihzh[IE][JE][ka];
static int ipos[LL]={0,0,0,0,0,0,0},jpos[LL]={2,10,15,40,50,60,70},kpos[LL]={0,0,0,0,0,0,0};
static int iwid,kwid,l,jgrid_end;
static float R_dist[IE][KE][LL],r_y[IE][KE][LL],N_shift[IE][KE][LL];
static float ezml[IE][KE],Eana[LL][TT],del_e[TT];
ic=IE/2;
jc=JE/2;
kc=KE/2;
ddx=.01; /* Set the cell size to 0.01 m */
dt=ddx/(2*3e8); /* Calculate the time step */
ib=IE-ia-10;
jb=JE-ja-10;
kb=KE-ka-10;
jmetal=50;
iwid=2;
kwid=4;
/*Initialize the arrays */
for(l=0;l<LL;l++){
for(i=0;i<=TT;i++){
Eana[l][i]=0.0;
}
}
for(j=0;j<JE;j++){
ez_inc[j]=0.0;
hx_inc[j]=0.0;
for(k=0;k<KE;k++){
for(i=0;i<IE;i++){
ex[i][j][k]=0.0;
ey[i][j][k]=0.0;
ez[i][j][k]=0.0;
dx[i][j][k]=0.0;
hz[i][j][k]=0.0;
gax[i][j][k]=1.0;
gay[i][j][k]=1.0;
gaz[i][j][k]=1.0;
}
}
}
for(i=0;i<ia;i++){
for(j=0;j<JE;j++){
for(k=0;k<KE;k++){
idxl[i][j][k]=0.0;
idxh[i][j][k]=0.0;
ihxl[i][j][k]=0.0;
ihxh[i][j][k]=0.0;
}
}
}
for(i=0;i<IE;i++){
for(j=0;j<ja;j++){
for(k=0;k<KE;k++){
idyl[i][j][k]=0.0;
idyh[i][j][k]=0.0;
ihyl[i][j][k]=0.0;
ihyh[i][j][k]=0.0;
}
}
}
for(i=0;i<IE;i++){
for(j=0;j<JE;j++){
for(k=0;k<ka;k++){
idzl[i][j][k]=0.0;
idzh[i][j][k]=0.0;
ihzl[i][j][k]=0.0;
ihzh[i][j][k]=0.0;
}
}
}
for(l=0;l<LL;l++){
for(i=ic-iwid;i<=ic+iwid;i++){
for(k=kc-kwid;k<=kc+kwid;k++){
R_dist[i][k][l]=sqrt(pow((ic+ipos[l]-i),2.0)
+pow((kc+kpos[l]-k),2.0)+pow((jpos[l]),2.0));
r_y[i][k][l]=jpos[l]/R_dist[i][k][l];
N_shift[i][k][l]=2*R_dist[i][k][l];
}
}
}
/* Bounday Conditions */
for(i=0;i<IE;i++){
gi1[i]=0.0;
gi2[i]=1.0;
gi3[i]=1.0;
fi1[i]=0.0;
fi2[i]=1.0;
fi3[i]=1.0;
}
for(j=0;j<JE;j++){
gj1[j]=0.0;
gj2[j]=1.0;
gj3[j]=1.0;
fj1[j]=0.0;
fj2[j]=1.0;
fj3[j]=1.0;
}
for(k=0;k<KE;k++){
gk1[k]=0.0;
gk2[k]=1.0;
gk3[k]=1.0;
fk1[k]=0.0;
fk2[k]=1.0;
fk3[k]=1.0;
}
printf("Number of PML cells--->");
scanf("%d",&npml);
printf("%d \n",npml);
n_pml=npml;
for(i=0;i<=n_pml;i++){
xxn =(npml-i)/npml;
xn =0.33*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
fi1[i]=xn;
fi1[IE-i-1]=xn;
gi2[i] =1.0/(1.0+xn);
gi2[IE-1-i] =1.0/(1.0+xn);
gi3[i] =(1.0-xn)/(1.0+xn);
gi3[IE-i-1] =(1.0-xn)/(1.0+xn);
xxn=(npml-i-0.5)/npml;
xn =0.33*pow(xxn,3.0);
gi1[i]=xn;
gi1[IE-i-2]=xn;
fi2[i] =1.0/(1.0+xn);
fi2[IE-2-i]=1.0/(1.0+xn);
fi3[i] =(1.0-xn)/(1.0+xn);
fi3[IE-2-i] =(1.0-xn)/(1.0+xn);
}
for(j=0;j<=n_pml;j++){
xxn =(npml-j)/npml;
xn =0.33*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
fj1[j]=xn;
fj1[JE-j-1]=xn;
gj2[j] =1.0/(1.0+xn);
gj2[JE-1-j] =1.0/(1.0+xn);
gj3[j] =(1.0-xn)/(1.0+xn);
gj3[JE-j-1] =(1.0-xn)/(1.0+xn);
xxn=(npml-j-0.5)/npml;
xn =0.33*pow(xxn,3.0);
gj1[j]=xn;
gj1[JE-j-2]=xn;
fj2[j] =1.0/(1.0+xn);
fj2[JE-2-j]=1.0/(1.0+xn);
fj3[j] =(1.0-xn)/(1.0+xn);
fj3[JE-2-j] =(1.0-xn)/(1.0+xn);
}
for(k=0;k<=n_pml;k++){
xxn =(npml-k)/npml;
xn =0.33*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
fk1[k]=xn;
fk1[KE-k-1]=xn;
gk2[k] =1.0/(1.0+xn);
gk2[KE-1-k] =1.0/(1.0+xn);
gk3[k] =(1.0-xn)/(1.0+xn);
gk3[KE-k-1] =(1.0-xn)/(1.0+xn);
xxn=(npml-k-0.5)/npml;
xn =0.33*pow(xxn,3.0);
gk1[k]=xn;
gk1[KE-k-2]=xn;
fk2[k] =1.0/(1.0+xn);
fk2[KE-2-k]=1.0/(1.0+xn);
fk3[k] =(1.0-xn)/(1.0+xn);
fk3[KE-2-k] =(1.0-xn)/(1.0+xn);
}
/* add a metal plane at j=jgrid_end*/
for(i=0;i<IE;i++){
for(k=0;k<KE;k++){
j=jmetal;
gax[i][j][k]=0.0;
gay[i][j][k]=0.0;
gaz[i][j][k]=0.0;
}
}
/* add a hole at ic,jc,kc*/
for(i=ic-iwid;i<=ic+iwid;i++){
for(k=kc-kwid;k<=kc+kwid;k++){
j=jmetal;
gax[i][j][k]=1.0;
gay[i][j][k]=1.0;
gaz[i][j][k]=1.0;
}
}
t0=40; /*Center of the incident pulse*/
spread=6.0; /*Width of the incident pulse*/
T=0;
NSTEPS=1;
fptime=fopen("time.dat","w");
fpt1=fopen("hou2.dat","w");
fpt2=fopen("hou10.dat","w");
fpt3=fopen("yuanchang15.dat","w");
fpt4=fopen("yuanchang40.dat","w");
fpt5=fopen("yuanchang50.dat","w");
fpt6=fopen("yuanchang60.dat","w");
fpt7=fopen("yuanchang70.dat","w");
while (NSTEPS>0) {
printf("NSTEPS-->"); /*NSTEPS is the number of times the*/
scanf("%d",&NSTEPS); /* main loop has executed*/
printf("%d\n",NSTEPS);
for (n=1; n<=NSTEPS; n++)
{
T=T+1; /*T keeps track of the total number*/
/* Start of the main FDTD Loop */
/* Calculate the incident buffer*/
for(j=1;j<JE;j++){
ez_inc[j]=ez_inc[j]+0.5*(hx_inc[j-1]-hx_inc[j]);
}
ez_inc[jmetal]=0.0;
for(i=ic-iwid; i<=ic+iwid; i++){
for(k=kc-kwid; k<=kc+kwid; k++){
ezml[i][k]=ez[i][jmetal][k];
}
}
/* Source */
pulse = sin(2*pi*1500*1e6*dt*T);
pulse = exp(-0.5*(pow((t0-T)/spread,2.0)));
ez_inc[3]=ez_inc[3]+pulse;
/* ABC for the incident buffer*/
ez_inc[0] =ez_low_m2;
ez_low_m2 =ez_low_m1;
ez_low_m1 =ez_inc[1];
ez_inc[JE-1] =ez_high_m2;
ez_high_m2 =ez_high_m1;
ez_high_m1 =ez_inc[JE-2];
/* Calculate the Dx field */
for(i=1;i<ia;i++){
for(j=jmetal+1;j<JE;j++){
for(k=1;k<KE;k++){
curl_h=(hx[i][j][k]-hx[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] = dx[i][j][k]+gj2[j]*gk2[k]*(curl_h+
gi1[i]*idxl[i][j][k]);
}
}
}
for(i=ia;i<=ib;i++){
for(j=jmetal+1;j<JE;j++){
for(k=1;k<KE;k++){
curl_h=(hz[i][j][k]-hz[i][j-1][k]-hy[i][j][k]+hy[i][j][k-1]);
dx[i][j][k] = gj3[j]dx[i][j][k]+gj2[j]*gk2[k]*0.5*curl_h;
}
}
}
/* Calculate 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]-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=jmetal+1;j<=jb;j++){
for(k=1;k<KE;k++){
curl_h=(hx[i][j][k]-hx[i][j][k-1]-hz[i][j][k]+hz[i-1][j][k]);
dy[i][j][k] = gi3[j]*gk3[j]*dy[i][j][j]+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]-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 Incident Dy */
for(i=ia;i<=ib;i++){
for(j=ja;j<=jb-1;j++){
dy[i][j][ka]=dy[i][j][ka]-0.5*hx_inc[j];
}
}
/* Calculate the Dz field */
for(i=1;i<IE;i++){
for(j=1;j<JE;j++){
for(k=0;k<ka;k++){
curl_h=(hy[i][j][k]-hy[i-1][j][k]-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++){
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -