一维fdtd模拟程序 C语言

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NUM_of_ZAXIS 400
float gauss_pulse(float T,float t0,float spread);
int file_save(float* data,char* filename);
void main()
{
float ex[NUM_of_ZAXIS],hy[NUM_of_ZAXIS];
float obj_parameters[NUM_of_ZAXIS][4]; /*模型参数设置*/
float CA[NUM_of_ZAXIS],CB[NUM_of_ZAXIS],CP[NUM_of_ZAXIS],CQ[NUM_of_ZAXIS]; float ca,cb,cp,cq; /*真空时的参数*/
float var_ca,var_cb,var_cp,var_cq; /* */
float
var_border,ex_low_m1,ex_low_m2,ex_high_m1,ex_high_m2,
ex_low_s1,ex_low_s2,ex_high_s1,ex_high_s2; /*边界吸收参数*/
float Epsilon,Mu,Pi,C; /*介电系数Epsilon 0，磁导系数Mu 0*/
float rel_epsz,rel_mu; /*相对介电系数，磁导系数*/
float e_sigma,h_sigma; /*电导率，磁导率*/
float dt,ddz;
float source,T;
int pos_driv_source;
int i,n,Nsteps;
int e_low,e_high;
FILE* fp;
/*初始化各个变量*/
Pi=3.14159;
Epsilon=8.85e-12;
Mu=4*Pi*1e-7;
C=pow( (float)(1/(Epsilon*Mu)),(float)(0.5));
ddz=0.01;
dt=ddz/(2*C);
T=0;
// 采用葛德彪书中一维FDTD公式中的参数e_sigma=0; h_sigma=0;
var_ca=0.5*e_sigma*dt/Epsilon;
var_cb=dt/Epsilon;
var_cp=0.5*h_sigma*dt/Mu;
var_cq=dt/Mu;
ca=(1-var_ca)/(1+var_ca);
cb=var_cb/(1+var_ca);
cp=(1-var_cp)/(1+var_cp);
cq=var_cq/(1+var_cp);
var_border=(C*dt-ddz)/(C*dt+ddz);
// 边界处用来存数据的临时变量
ex_low_m1=0;
ex_low_m2=0;
ex_high_m1=0;
ex_high_m2=0;
ex_low_s1=0;
ex_low_s2=0;
ex_high_s1=0;
ex_high_s2=0;
// 上，下边界位置
e_low=20;
e_high=380;
pos_driv_source=200;
for(i=0;i<NUM_of_ZAXIS;i++) {
ex=0;
hy=0;
obj_parameters[0]=1;
obj_parameters[1]=1;
obj_parameters[2]=0;
obj_parameters[3]=0; CA=ca;
CB=cb;
CP=cp;
CQ=cq;
}
// 开始主循环
Nsteps=1;
while(Nsteps>0)
{
printf("Nsteps-->");
scanf("%d",&Nsteps);
for(n=0;n<Nsteps;n++)
{
T+=1;
// 计算Ex
for(i=e_low+1;i<e_high;i++)
{
// 采用葛德彪书中的一维FDTD 公式ex=CA*ex-CB*(hy-hy[i-1])/ddz; }
if(T<200){
ex[pos_driv_source]+=gauss_pulse(T,40,12);
}
//Mur 边界吸收
ex_low_m2=ex_low_m1;
ex_low_m1=ex[e_low];
ex_low_s2=ex_low_s1;
ex_low_s1=ex[e_low+1];
ex_high_m2=ex_high_m1;
ex_high_m1=ex[e_high];
ex_high_s2=ex_high_s1;
ex_high_s1=ex[e_high-1];
ex[e_low]=ex_low_s2+var_border*(ex[e_low+1]-ex_low_m2); ex[e_high]=ex_high_s2+var_border*(ex[e_high-1]-ex_high_m2);
// 行波时延法
// ex[e_high]=0.5*(ex_high_m2+ex_high_s2);
// ex[e_low]=0.5*(ex_low_m2+ex_low_s2);
// 计算Hy
for(i=e_low;i<e_high;i++)
{
hy=CP*hy-CQ*(ex[i+1]-ex)/ddz;
}
}
// 存储数据
fp=fopen("Ex","w");
for(i=0;i<NUM_of_ZAXIS;i++)
{
fprintf(fp,"%6.2f \n",ex);
}
fclose(fp);
fp=fopen("Hy","w");
for(i=0;i<NUM_of_ZAXIS;i++)
{
fprintf(fp,"%6.5f \n",hy);
}
fclose(fp);
printf("T=%6.2f",T);
}
}
float gauss_pulse(float T,float t0,float spread) {
float pulse;
pulse=exp(-0.5*(pow((T-t0)/spread,2))); return pulse;
}。