几个用Matlab和C写的关于FDTD程序:希望对大家能有所帮助！
//fd2d_3.1.c 2D TM program
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define IE 60
#define JE 60
main()
{
float ga[IE][JE], dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE];
int l, n, i, j, ic, jc, nsteps;
float ddx, dt, T, epsz, pi, epsilon, sigma, eaf;
float t0, spread, pulse;
FILE *fp;
ic=IE/2;
jc=JE/2;
ddx=0.01;
dt=ddx/6e8;
epsz=8.8e-12;
pi=3.14159;
for (j=0; j<JE; j++){
printf("%2d ",j);
for (i=0;i<IE;i++){
dz[j]=0.0;
ez[j]=0.0;
hx[j]=0.0;
hy[j]=0.0;
ga[j]=1.0;
printf("%5.2f ", ga[j]);
}
printf("\n");
}
t0=30;
spread=6.0;
T=0;
nsteps=1;
while (nsteps>0){
printf("nsteps --> ");
scanf("%d", &nsteps);
for (n=1; n<=nsteps; n++){
T=T+1;
/* start of the Main loop */
//calculate the Dz field
for(j=1;j<JE;j++){
for(i=1;i<IE;i++){
dz[j]=dz[j]+0.5*(hy[j]-hy[i-1][j]-hx[j]+hx[j-1]);
}
}
//put a Guassian pulse in thye middle
pulse=exp(-0.5*pow( (t0-T)/spread,2.0 ) );
dz[ic][jc]=pulse;
//calculate the Ez field
for(j=0;j<JE;j++){
for(i=0;i<IE;i++){
ez[j]=ga[j]*dz[j];
}
}
//calculate the Hx field
for(j=0; j<JE-1;j++){
for(i=0;i<IE;i++){
hx[j]=hx[j]+0.5*(ez[j]-ez[j+1]);
}
}
//calculate the Hy field
for(j=0; j<JE;j++){
for(i=0;i<IE-1;i++){
hy[j]=hy[j]+0.5*(ez[i+1][j]-ez[j]);
}
}
}
/*-----------End of the main FDTD loop-------------*/
for(j=1;j<jc;j++){
printf("%2d ", j);
for(i=1;i<ic;i++){
printf("%5.2f ",ez[2*i][2*j]);
}
printf("\n");
}
//printf("T=%5.0f \n",T);
/* write the E fireld to a file "Ez.xls" */
fp=fopen("Ez.xls", "w");
for(j=0;j<JE;j++){
for(i=0;i<IE;i++){
fprintf(fp, "%6.3f \t",ez[j]);
}
fprintf(fp,"\n");
}
fclose(fp);
printf("T=%5.0f \n",T);
}
}
//fd2d_3.2.c 2D TM program
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define IE 60
#define JE 60
main()
{
float ga[IE][JE], dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE];
int l, n, i, j, ic, jc, nsteps, npml;
float ddx, dt, T, epsz, pi, epsilon, sigma, eaf;
float xn, xxn, xnum, xd, curl_e;
float t0, spread, pulse;
float gi2[IE], gi3[IE];
float gj2[JE], gj3[JE];
float fi1[IE],fi2[IE],fi3[IE];
float fj1[JE],fj2[JE],fj3[JE];
float ihx[IE][JE], ihy[IE][JE];
FILE *fp;
ic=IE/2-15;
jc=JE/2-15;
ddx=0.01;
dt=ddx/6e8;
epsz=8.8e-12;
pi=3.14159;
for (j=0; j<JE; j++){
printf("%2d ",j);
for (i=0;i<IE;i++){
dz[j]=0.0;
ez[j]=0.0;
hx[j]=0.0;
hy[j]=0.0;
ihx[j]=0.0;
ihy[j]=0.0;
ga[j]=1.0;
printf("%5.2f ", ga[j]);
}
printf("\n");
}
/*----------calculate the PML parameters-------------------*/
for(i=0;i<IE;i++){
gi2=1.0;
gi3=1.0;
fi1=0.0;
fi2=1.0;
fi3=1.0;
}
for(j=0;j<JE;j++){
gj2[j]=1.0;
gj3[j]=1.0;
fj1[j]=0.0;
fj2[j]=1.0;
fj3[j]=1.0;
}
printf("Number of PML cells --> ");
scanf("%d", &npml);
for(i=0;i<=npml;i++){
xnum=npml-i;
xd=npml;
xxn=xnum/xd;
xn=0.33*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n", i, xxn, xn);
gi2=1.0/(1.0+xn);
gi2[IE-1-i]=1.0/(1.0+xn);
gi3=(1.0-xn)/(1.0+xn);
gi3[IE-1-i]=(1.0-xn)/(1.0+xn);
xxn=(xnum-0.5)/xd; //orig
//xxn=(xnum)/xd;
xn=0.25*pow(xxn,3.0);
fi1=xn;
fi1[IE-2-i]=xn;
fi2=1.0/(1.0+xn);
fi2[IE-2-i]=1.0/(1.0+xn);
fi3=(1.0-xn)/(1.0+xn);
fi3[IE-2-i]=(1.0-xn)/(1.0+xn);
}
for(j=0;j<=npml;j++){
xnum=npml-j;
xd=npml;
xxn=xnum/xd;
xn=0.33*pow(xxn,3.0);
printf("%d %7.4f %7.4f \n", i, xxn, 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-1-j]=(1.0-xn)/(1.0+xn);
xxn=(xnum-0.5)/xd; //orig
//xxn=(xnum)/xd;
xn=0.25*pow(xxn,3.0);
fj1[j]=xn;
fj1[JE-2-j]=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);
}
t0=40;
spread=15.0;
T=0;
nsteps=1;
while (nsteps>0){
printf("nsteps --> ");
scanf("%d", &nsteps);
for (n=1; n<=nsteps; n++){
T=T+1;
/* start of the Main loop */
//calculate the Dz field
for(j=1;j<JE;j++){
for(i=1;i<IE;i++){
dz[j]=gi3*gj3[j]*dz[j]
+gi2*gj2[j]*0.5*(hy[j]-hy[i-1][j]-hx[j]+hx[j-1]);
}
}
//put a Guassian pulse in thye middle
//pulse=exp(-0.5*pow( (t0-T)/spread,2.0 ) );
pulse=sin(2*pi*1500*1e6*dt*T);
dz[ic][jc]=pulse;
//dz[ic][jc]=dz[ic][jc]+pulse;
//(1)
//calculate the Ez field
for(j=1;j<JE-1;j++){
//for(i=0;i<IE;i++){
for(i=1;i<IE-1;i++){
ez[j]=ga[j]*dz[j];
}
}
//(2)
/* set the Ez edges to 0, as part of the PML */
for(j=0;j<JE;j++){
ez[0][j]=0.0;
ez[IE-1][j]=0.0;
}
for(i=0;i<IE;i++){
ez[0]=0.0;
ez[JE-1]=0.0;
}
//(3)
//calculate the Hx field
for(j=0; j<JE;j++){
for(i=0;i<IE;i++){
curl_e=ez[j]-ez[j+1];
ihx[j]=ihx[j]+fi1*curl_e;
hx[j]=fj3[j]*hx[j]
+fj2[j]*0.5*(curl_e+ihx[j]);
}
}
//calculate the Hy field
for(j=0; j<JE;j++){
for(i=0;i<IE;i++){
curl_e=ez[i+1][j]-ez[j];
ihy[j]=ihy[j]+fj1[j]*curl_e;
hy[j]=fi3*hy[j]
+fi2*0.5*(curl_e+ihy[j]);
}
}
}
/*-----------End of the main FDTD loop-------------*/
for(j=1;j<jc;j++){
printf("%2d ", j);
for(i=1;i<ic;i++){
printf("%5.2f ",ez[2*i][2*j]);
}
printf("\n");
}
//printf("T=%5.0f \n",T);
/* write the E fireld to a file "Ez.xls" */
fp=fopen("Ez.xls", "w");
for(j=0;j<JE;j++){
for(i=0;i<IE;i++){
fprintf(fp, "%6.3f \t",ez[j]);
}
fprintf(fp,"\n");
}
fclose(fp);
printf("T=%5.0f \n",T);
}
}
//fd2d_3.3.c 2D TM program with plane wave source
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define IE 60
#define JE 60
main()
{
float ez_inc[JE], hx_inc[JE];
float ez_inc_low_m1, ez_inc_low_m2;
float ez_inc_high_m1, ez_inc_high_m2;
int ia, ib, ja, jb;
float ga[IE][JE], dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE];
int l, n, i, j, ic, jc, nsteps, npml;
float ddx, dt, T, epsz, pi, epsilon, sigma, eaf;
float xn, xxn, xnum, xd, curl_e;
float t0, spread, pulse;
float gi2[IE], gi3[IE];
float gj2[JE], gj3[JE];
float fi1[IE],fi2[IE],fi3[IE];
float fj1[JE],fj2[JE],fj3[JE];
float ihx[IE][JE], ihy[IE][JE];
FILE *fp;
ic=IE/2-15;
jc=JE/2-15;
ia=7; //Total/scattered field boundaries
ib=IE-ia-1;
ja=7;
jb=JE-ja-1;
ddx=0.01; //cell size
dt=ddx/6e8; //time step
epsz=8.8e-12;
pi=3.14159;
for (j=0; j<JE; j++){
printf("%2d ",j);
for (i=0;i<IE;i++){
dz[j]=0.0;
ez[j]=0.0;
hx[j]=0.0;
hy[j]=0.0;
ihx[j]=0.0;
ihy[j]=0.0;
ga[j]=1.0;
printf("%5.2f ", ga[j]);
}
printf("\n");
}
for (j=0; j<JE; j++){
ez_inc[j]=0, hx_inc[j]=0;
}
ez_inc_low_m1=0, ez_inc_low_m2=0;
ez_inc_high_m1=0, ez_inc_high_m2=0;
/*----------calculate the PML parameters-------------------*/
for(i=0;i<IE;i++){
gi2=1.0;
gi3=1.0;
fi1=0.0;
fi2=1.0;
fi3=1.0;
}
for(j=0;j<JE;j++){
gj2[j]=1.0;
gj3[j]=1.0;
fj1[j]=0.0;
fj2[j]=1.0;
fj3[j]=1.0;
}
printf("Number of PML cells --> ");
scanf("%d", &npml);
for(i=0;i<=npml;i++){
xnum=npml-i;
xd=npml;
xxn=xnum/xd;
xn=0.33*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n", i, xxn, xn);
gi2=1.0/(1.0+xn);
gi2[IE-1-i]=1.0/(1.0+xn);
gi3=(1.0-xn)/(1.0+xn);
gi3[IE-1-i]=(1.0-xn)/(1.0+xn);
xxn=(xnum-0.5)/xd; //orig
//xxn=(xnum)/xd;
xn=0.25*pow(xxn,3.0);
fi1=xn;
fi1[IE-2-i]=xn;
fi2=1.0/(1.0+xn);
fi2[IE-2-i]=1.0/(1.0+xn);
fi3=(1.0-xn)/(1.0+xn);
fi3[IE-2-i]=(1.0-xn)/(1.0+xn);
}
for(j=0;j<=npml;j++){
xnum=npml-j;
xd=npml;
xxn=xnum/xd;
xn=0.33*pow(xxn,3.0);
printf("%d %7.4f %7.4f \n", j, xxn, 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-1-j]=(1.0-xn)/(1.0+xn);
xxn=(xnum-0.5)/xd; //orig
//xxn=(xnum)/xd;
xn=0.25*pow(xxn,3.0);
fj1[j]=xn;
fj1[JE-2-j]=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);
}
t0=20;
spread=8.0;
T=0;
nsteps=1;
while (nsteps>0){
printf("nsteps --> ");
scanf("%d", &nsteps);
for (n=1; n<=nsteps; n++){
T=T+1;
/* start of the Main loop */
for (j=1; j<JE-1;j++){
ez_inc[j]=ez_inc[j]+0.5*(hx_inc[j-1]-hx_inc[j]);
}
//ABC for the incident buffer
ez_inc[0] =ez_inc_low_m2;
ez_inc_low_m2=ez_inc_low_m1;
ez_inc_low_m1=ez_inc[1];
//ez_inc_low_m2=ez_inc_low_m1;
//ez_inc[0] =ez_inc_low_m2;
ez_inc[JE-1] =ez_inc_high_m2;
ez_inc_high_m2=ez_inc_high_m1;
ez_inc_high_m1=ez_inc[JE-2];
//ez_inc_high_m2=ez_inc_high_m1;
//ez_inc[JE-1] =ez_inc_high_m2;
//calculate the Dz field
for(j=1;j<JE;j++){
for(i=1;i<IE;i++){
dz[j]=gi3*gj3[j]*dz[j]
+gi2*gj2[j]*0.5*(hy[j]-hy[i-1][j]
-hx[j]+hx[j-1]);
}
}
//source
pulse=exp(-0.5*( pow((t0-T)/spread,2.0) ));
ez_inc[3]=pulse;
//Incident Dz value
for (i=ia;i<=ib;i++){
dz[ja]=dz[ja]+0.5*hx_inc[ja-1];
dz[jb]=dz[jb]-0.5*hx_inc[jb];
}
//(1)
//calculate the Ez field
for(j=0;j<JE;j++){
//for(i=0;i<IE;i++){
for(i=0;i<IE;i++){
ez[j]=ga[j]*dz[j];
}
}
//(2)
/* set the Ez edges to 0, as part of the PML */
for(j=0;j<JE-1;j++){
ez[0][j]=0.0;
ez[IE-1][j]=0.0;
}
for(i=0;i<IE-1;i++){
ez[0]=0.0;
ez[JE-1]=0.0;
}
for(j=0;j<JE;j++){
hx_inc[j]=hx_inc[j]+0.5*(ez_inc[j]-ez_inc[j+1]);
}
//(3)
//calculate the Hx field
for(j=0; j<JE-1;j++){
for(i=0;i<IE;i++){
curl_e=ez[j]-ez[j+1];
ihx[j]=ihx[j]+fi1*curl_e;
hx[j]=fj3[j]*hx[j]
+fj2[j]*0.5*(curl_e+ihx[j]);
}
}
//incident Hx value
for (i=ia;i<=ib;i++){
hx[ja-1] =hx[ja-1] + 0.5*ez_inc[ja];
hx[jb] =hx[jb] - 0.5*ez_inc[jb];
}
//calculate the Hy field
for(j=0; j<JE;j++){
for(i=0;i<IE;i++){
curl_e=ez[i+1][j]-ez[j];
ihy[j]=ihy[j]+fj1[j]*curl_e;
hy[j]=fi3*hy[j]
+fi2*0.5*(curl_e+ihy[j]);
}
}
//incident Hy value
for (j=ja;j<=jb;j++){
hy[ia-1][j] =hy[ia-1][j] - 0.5*ez_inc[j];
hy[ib][j] =hy[ib][j] + 0.5*ez_inc[j];
}
}
/*-----------End of the main FDTD loop-------------*/
for(j=1;j<jc;j++){
printf("%2d ", j);
for(i=1;i<ic;i++){
printf("%5.2f ",ez[2*i][2*j]);
}
printf("\n");
}
//printf("T=%5.0f \n",T);
/* write the E fireld to a file "Ez.xls" */
fp=fopen("Ez.xls", "w");
for(j=0;j<JE;j++){
for(i=0;i<IE;i++){
fprintf(fp, "%6.3f \t",ez[j]);
}
fprintf(fp,"\n");
}
fclose(fp);
printf("T=%5.0f \n",T);
}
}
//fd2d_3.4.c 2D TM simulation of a plane
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define IE 60
#define JE 60
#define NFREQS 3
main()
{
float real_pt[NFREQS][IE][JE], imag_pt[NFREQS][IE][JE];
float real_amp[IE][JE], imag_amp[IE][JE];
float real_in[5], imag_in[5], amp_in[5], phase_in[5], freq[3], arg[3];
float ez_inc[JE], hx_inc[JE];
float ez_inc_low_m1, ez_inc_low_m2;
float ez_inc_high_m1, ez_inc_high_m2;
int ia, ib, ja, jb;
float ga[IE][JE], gb[IE][JE], dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE], iz[IE][JE];
int l, n, m, i, j, ic, jc, nsteps, npml;
float ddx, dt, T, epsz, pi, epsilon, sigma, eaf;
float xn, xxn, xnum, xd, curl_e;
float t0, spread, pulse;
float gi2[IE], gi3[IE];
float gj2[JE], gj3[JE];
float fi1[IE],fi2[IE],fi3[IE];
float fj1[JE],fj2[JE],fj3[JE];
float ihx[IE][JE], ihy[IE][JE];
FILE *fp;
ic=IE/2;
jc=JE/2;
ia=7; //Total/scattered field boundaries
ib=IE-ia-1;
ja=7;
jb=JE-ja-1;
ddx=0.01; //cell size
dt=ddx/6e8; //time step
epsz=8.8e-12;
pi=3.14159;
for (j=0; j<JE; j++){
printf("%2d ",j);
for (i=0;i<IE;i++){
dz[j]=0.0;
ez[j]=0.0;
hx[j]=0.0;
hy[j]=0.0;
ihx[j]=0.0;
ihy[j]=0.0;
iz[j]=0.0;
ga[j]=1.0;
gb[j]=0.0;
printf("%5.2f ", ga[j]);
}
printf("\n");
}
for (j=0; j<JE; j++){
ez_inc[j]=0, hx_inc[j]=0;
}
ez_inc_low_m1=0, ez_inc_low_m2=0;
ez_inc_high_m1=0, ez_inc_high_m2=0;
/*----------calculate the PML parameters-------------------*/
for(i=0;i<IE;i++){
gi2=1.0;
gi3=1.0;
fi1=0.0;
fi2=1.0;
fi3=1.0;
}
for(j=0;j<JE;j++){
gj2[j]=1.0;
gj3[j]=1.0;
fj1[j]=0.0;
fj2[j]=1.0;
fj3[j]=1.0;
}
printf("Number of PML cells --> ");
scanf("%d", &npml);
for(i=0;i<=npml;i++){
xnum=npml-i;
xd=npml;
xxn=xnum/xd;
xn=0.33*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n", i, xxn, xn);
gi2=1.0/(1.0+xn);
gi2[IE-1-i]=1.0/(1.0+xn);
gi3=(1.0-xn)/(1.0+xn);
gi3[IE-1-i]=(1.0-xn)/(1.0+xn);
xxn=(xnum-0.5)/xd; //orig
//xxn=(xnum)/xd;
xn=0.25*pow(xxn,3.0);
fi1=xn;
fi1[IE-2-i]=xn;
fi2=1.0/(1.0+xn);
fi2[IE-2-i]=1.0/(1.0+xn);
fi3=(1.0-xn)/(1.0+xn);
fi3[IE-2-i]=(1.0-xn)/(1.0+xn);
}
for(j=0;j<=npml;j++){
xnum=npml-j;
xd=npml;
xxn=xnum/xd;
xn=0.33*pow(xxn,3.0);
printf("%d %7.4f %7.4f \n", j, xxn, 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-1-j]=(1.0-xn)/(1.0+xn);
xxn=(xnum-0.5)/xd; //orig
//xxn=(xnum)/xd;
xn=0.25*pow(xxn,3.0);
fj1[j]=xn;
fj1[JE-2-j]=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);
}
/*parameters for Fourier transform*/
freq[0]=50e6;
freq[1]=300e6;
freq[2]=700e6;
for(n=0;n<NFREQS;n++){
arg[n]=2*pi*freq[n]*dt;
}
/*specify the dielectric cylinder*/
float radius, dist, xdist, ydist;
printf("Cylinder radius (cells), epsilon, sigma --->");
scanf("%f %f %f", &radius, &epsilon, &sigma);
printf("Radius=%f, Eps=%f, sigma=%f \n",radius, epsilon, sigma);
for(j=ja;j<jb;j++){
for(i=ia;i<ib;i++){
xdist=(ic-i);
ydist=(jc-j);
dist=sqrt( pow(xdist,2.0)+pow(ydist,2.0) );
if(dist<=radius){
ga[j]=1.0/(epsilon+(sigma*dt/epsz));
gb[j]=sigma*dt/epsz;
}
}
}
printf(" Ga");
for(j=ja;j<jb;j++){
for(i=ia;i<ib;i++){
printf("%5.2f",ga[j]);
}
printf("\n");
}
printf(" Gb");
for(j=ja;j<jb;j++){
for(i=ia;i<ib;i++){
printf("%5.2f",gb[j]);
}
printf("\n");
}
t0=20;
spread=8.0;
T=0;
nsteps=1;
while (nsteps>0){
printf("nsteps --> ");
scanf("%d", &nsteps);
for (n=1; n<=nsteps; n++){
T=T+1;
/* start of the Main loop */
//calculate the incident Ez
for (j=1; j<JE;j++){
ez_inc[j]=ez_inc[j]+0.5*(hx_inc[j-1]-hx_inc[j]);
}
/*Fourier transform of the incident field*/
for(m=0;m<NFREQS;m++){
real_in[m]=real_in[m]+cos(arg[m]*T)*ez_inc[ja-1];
imag_in[m]=imag_in[m]-sin(arg[m]*T)*ez_inc[ja-1];
}
//ABC for the incident buffer
ez_inc[0] =ez_inc_low_m2;
ez_inc_low_m2=ez_inc_low_m1;
ez_inc_low_m1=ez_inc[1];
//ez_inc_low_m2=ez_inc_low_m1;
//ez_inc[0] =ez_inc_low_m2;
ez_inc[JE-1] =ez_inc_high_m2;
ez_inc_high_m2=ez_inc_high_m1;
ez_inc_high_m1=ez_inc[JE-2];
//ez_inc_high_m2=ez_inc_high_m1;
//ez_inc[JE-1] =ez_inc_high_m2;
//calculate the Dz field
for(j=1;j<JE;j++){
for(i=1;i<IE;i++){
dz[j]=gi3*gj3[j]*dz[j]
+gi2*gj2[j]*0.5*(hy[j]-hy[i-1][j]
-hx[j]+hx[j-1]);
}
}
//source
pulse=exp(-0.5*( pow((t0-T)/spread,2.0) ));
ez_inc[3]=pulse;
//Incident Dz value
for (i=ia;i<=ib;i++){
dz[ja]=dz[ja]+0.5*hx_inc[ja-1];
dz[jb]=dz[jb]-0.5*hx_inc[jb];
}
//(1)
//calculate the Ez field
for(j=0;j<JE;j++){
//for(i=0;i<IE;i++){
for(i=0;i<IE;i++){
ez[j]=ga[j]*(dz[j]-iz[j]);
iz[j]=iz[j]+gb[j]*ez[j];
}
}
/* calculate the fourier transform of Ex*/
for(j=0;j<JE;j++){
for(i=0;i<IE;i++){
for(m=0;m<NFREQS;m++){
real_pt[m][j]=real_pt[m][j]
+cos(arg[m]*T)*ez[j];
imag_pt[m][j]=imag_pt[m][j]
+sin(arg[m]*T)*ez[j];
}
}
}
//(2)
/* set the Ez edges to 0, as part of the PML */
for(j=0;j<JE-1;j++){
ez[0][j]=0.0;
ez[IE-1][j]=0.0;
}
for(i=0;i<IE-1;i++){
ez[0]=0.0;
ez[JE-1]=0.0;
}
//calculate the incident Hx
for(j=0;j<JE;j++){
hx_inc[j]=hx_inc[j]+0.5*(ez_inc[j]-ez_inc[j+1]);
}
//(3)
//calculate the Hx field
for(j=0; j<JE-1;j++){
for(i=0;i<IE;i++){
curl_e=ez[j]-ez[j+1];
ihx[j]=ihx[j]+fi1*curl_e;
hx[j]=fj3[j]*hx[j]
+fj2[j]*0.5*(curl_e+ihx[j]);
}
}
//incident Hx value
for (i=ia;i<=ib;i++){
hx[ja-1] =hx[ja-1] + 0.5*ez_inc[ja];
hx[jb] =hx[jb] - 0.5*ez_inc[jb];
}
//calculate the Hy field
for(j=0; j<JE;j++){
for(i=0;i<IE;i++){
curl_e=ez[i+1][j]-ez[j];
ihy[j]=ihy[j]+fj1[j]*curl_e;
hy[j]=fi3*hy[j]
+fi2*0.5*(curl_e+ihy[j]);
}
}
//incident Hy value
for (j=ja;j<=jb;j++){
hy[ia-1][j] =hy[ia-1][j] - 0.5*ez_inc[j];
hy[ib][j] =hy[ib][j] + 0.5*ez_inc[j];
}
}
/*-----------End of the main FDTD loop-------------*/
/*calculate the Fourier amplitude and phase of the incident pulse*/
for(j=1;j<jc;j++){
printf("%2d ", j);
for(i=1;i<ic;i++){
printf("%5.2f ",ez[2*i][2*j]);
}
printf("\n");
}
//printf("T=%5.0f \n",T);
/* write the E fireld to a file "Ez.xls" */
fp=fopen("Ez.xls", "w");
for(j=0;j<JE;j++){
for(i=0;i<IE;i++){
fprintf(fp, "%6.3f \t",ez[j]);
}
fprintf(fp,"\n");
}
fclose(fp);
printf("T=%5.0f \n",T);
}
}
%======================================================================================
% function Ex_obs = FDTD_1D(Ex,dx,dt,srcIndex,srcJx,pf,nSteps,e_r,mu_r,sigma,sigma_m,k_obs)
%
% This function performs one-dimensional (Cartesian) finite-difference time-domain
% simulations. It assumes the spatial grid extends from k = 1 to k = Nz and that the
% simulation runs for 'nSteps' time steps.
% The function should take the following input arguments, in this order:
%
% Input Parameter Size Description
% Ex 1xNz initial condition for Ex (V/m)
% dx 1x1 spatial step size (m)
% dt 1x1 time step size (s)
% Src_indx 1x1 grid index for source location
% Src_Jx 1xnSteps source current stimulus (A/m2)
% pf 1x1 plot flag - # steps between plots
% nSteps 1x1 # steps for simulation
% e_r 1xNz relative permittivity
% mu_r 1x(Nz-1) relative permeability
% sigma 1xNz conductivity (S/m)
% sigma_m 1x(Nz-1) magnetic conductivity (Wb/C穖)
% k_obs 1x1 index of observation point
% The one output parameter should be Ex_obs, a 1xnSteps vector
% of values at the observation point. If no observation point is desired,
% set k_obs = [ ] (an empty matrix) and return [ ] in Ex_obs.
% The function also generates a plot of Ex every 'pf' time steps
% so that we can observe what is going on
%
% Written by
% Aroh Barjatya
% For Utah State University class ECE 5830 Electromagnetics 2
function Ex_obs = FDTD_1D(Ex,dx,dt,srcIndex,srcJx,pf,nSteps,e_r,mu_r,sigma,sigma_m,k_obs)
c1 = ((2 .* e_r) - (sigma .* dt)) ./ ((2 .* e_r) + (sigma .* dt));
c2 = (2 .* dt) ./ (dx .* ((2 .* e_r) + (sigma .* dt)));
c3 = (2 .* dt) ./ ((2 .* e_r) + (sigma .* dt));
c4 = ((2 .* mu_r) - (sigma_m .* dt)) ./ ((2 .* mu_r) + (sigma_m .* dt));
c5 = (2 .* dt) ./ (dx .* ((2 .* mu_r) + (sigma_m .* dt)));
Ex(srcIndex) = srcJx(1);
Hy = zeros(1,length(Ex) - 1);
Ex_obs = zeros(1,nSteps);
if length(k_obs) > 0
Ex_obs(1) = Ex(k_obs);
end
for i = 2:nSteps
Ex(2:end-1) = c1(2:end-1) .* Ex(2:end-1) - c2(2:end-1) .* (Hy(2:end)-Hy(1:end-1));
Ex(srcIndex) = Ex(srcIndex) - c3(srcIndex)*srcJx(i);
Hy = (c4 .* Hy) - (c5 .* (Ex(2:end) - Ex(1:end-1)));
if (mod(i, pf) == 0)
plot(Ex);
a = sprintf('step number = %d', i);
title(a);
set(gcf,'color','white');
xlabel('Grid Points -->');
ylabel('Magnitude -->');
% axis([1 1001 -1 1]);
pause(0.2);
end
if length(k_obs) > 0
Ex_obs(i) = Ex(k_obs) ;
end
end
% function p=Gaussian_pulse(dt,f_c,len)
%
% The function creates a Gaussian pulse centered at frequency fc
% It assumes a sample frequency of fs = 1/dt
%
% the pulse vector will be of length len; specify too short of a value for len
% and an error will be generated; delay is needed before the start of the pulse
% to reduce the step function transient created at t=0
%
% the resulting pulse has a lowpass characteristic with a 3 dB bandwidth f_c
% (i.e. extending from roughly DC to f_c)
%
% Originally Written by Michael Tompkins, Asst. Prof at USU Fall 2004
% Minor modifications by Aroh Barjatya
function p=Gaussian_pulse(dt,f_c,len)
if 0
Gaussian_pulse(1e-7,2e5,15000);
Gaussian_pulse(1.6e-10,30e6,15000);
end
spread = 1/2/pi/f_c/dt; % this is the standard deviation of our Gaussian function
% must wait a while to launch the pulse so that we don't get too big of a
% "jump" (step function) at t=0; this step has higher frequency harmonics
% in it that will not meet our requirement that dx << lambda
delay = 5*spread;
if delay>len
error('Need longer simulation to accommodate desired pulse');
end
n = 0:len-1;
p = exp(-((n-delay)/spread).^2); % Gaussian pulse
%=====================================================================================
% This is a test script for FDTD_1D
% We test it by giving it a Gaussian-pulse
% plane-wave source in the middle of a 1001-point grid.
% Assumed free space everywhere.
%
% Written by
% Aroh Barjatya
% For Utah State University class ECE 5830 Electromagnetics 2
clear all
close all
% Physical Constants
e_0 = 8.85e-12;
mu_0 = 4e-7 * pi;
c = 2.99e8;
eta_0 = sqrt(mu_0/e_0);
% Simulation Constants
Nz = 1001;
nSteps = 3500;
Ex = zeros(1,Nz);
dx = 0.1;
dt = dx / (2*c);
e_r = ones(1,Nz) .* e_0;
mu_r = ones(1,Nz - 1) .* mu_0;
pf = 20;
% Define the source and its location
pulse = Gaussian_pulse(dt,30e6,nSteps);
srcJx = (-2/(eta_0*dx))*pulse;
srcIndex = (Nz - 1)/2;
sigma = zeros(1,Nz);
sigma_m = zeros(1,Nz-1);
FDTD_1D(Ex, dx, dt, srcIndex, srcJx, pf, nSteps, e_r, mu_r, sigma, sigma_m, []);
%====================================================================================
% This is the second test script for FDTD_1D
% In addition to what was done in test cript 1, we add a 50 unit think perfectly
% matched layer (PML) on both ends of the simulation space.
%
% Written by
% Aroh Barjatya
% For Utah State University class ECE 5830 Electromagnetics 2
clear all
close all
% Physical Constants
e_0 = 8.85e-12;
mu_0 = 4e-7 * pi;
c = 2.99e8;
eta_0 = sqrt(mu_0/e_0);
% Simulation Constants
Nz = 1001;
nSteps = 3500;
Ex = zeros(1,Nz);
dx = 0.1;
dt = dx / (2*c);
e_r = ones(1,Nz) .* e_0;
mu_r = ones(1,Nz - 1) .* mu_0;
pf = 20;
% Define the source and its location
pulse = Gaussian_pulse(dt,30e6,nSteps);
srcJx = (-2/(eta_0*dx))*pulse;
srcIndex = (Nz - 1)/2;
sigma = zeros(1,Nz);
sigma_m = zeros(1,Nz-1);
% for PML the following lines define sigma and sigma_m
% if no PML is required, then comment the following lines
% The Left Edge
k = 50;
k0 = 50;
for i = k:-1:1
sigma(i) = ((2 * e_r(i)) / dt) * 0.33 * ((k0 - i + 1) / k)^3;
sigma_m(i) = sigma(i) * mu_r(i) / e_r(i);
end
% The Right Edge
k0 = Nz-k;
for i = Nz-k:Nz
sigma(i) = ((2 * e_r(i)) / dt) * 0.33 * ((i - k0 + 1) / k)^3;
sigma_m(i-1) = sigma(i) * mu_r(i-1) / e_r(i);
end
FDTD_1D(Ex, dx, dt, srcIndex, srcJx, pf, nSteps, e_r, mu_r, sigma, sigma_m, []);
% 简化版，介质（非磁性）不能接触边界
N=6; % 运算步数
I=200; % 空间大小
J=400; % 空间大小
c=3e8; % 真空光速
L=632.8e-9;
w=2*pi*c/L; % 角频率=2.9788e+015
e0=(1e-9)/(36*pi); % 真空的介电常数=8.8419e-012
e=ones(I,J)*e0;
u0=pi*4e-7; % 真空的磁导率=1.2566e-006
ds=4e-9; % 波长=2*pi*c/w 步长=波长/40
dt=ds/(2*c); % 时间步
g0=0;
g=zeros(I,J); % 真空的电导率
gm=0; % 真空的导磁率
sourcePosX=50;
sourcePosY=40;
ArrayPosX=10;
ArrayPosY=100;
ArrayLength=100;
ArrayWidth=3;
period=6;
ArrayNumber=floor( (I-2*ArrayPosX)/period );
% 电磁场分量定义
Hz=zeros(I,J);
Ex=zeros(I,J);
Ey=zeros(I,J);
Hz1=zeros(I,J);
Ex1=zeros(I,J);
Ey1=zeros(I,J);
% Source -----------------------------------------------------------
Source=zeros(I,J);
%Source(sourcePosX,sourcePosY)=1;
x1=sourcePosX;
x2=I-sourcePosX*2;
for x=1:1:x2-x1+1;
Source(x+x1-1,sourcePosY)=exp(-4*pi*(x-(x2-x1)/2)^2/(x2-x1)^2);
end
% --------------------------------------------------------------------
% Object --------------------------------------------------------------
nn=L*1e6; % L is wave length
coon=(-6.21198035452752*nn.^5+ 36.9470450716619*nn.^4-83.1702487377279*nn.^3+...
94.5079035255191*nn.^2-48.61825472792*nn+9.64235361750187)
eee=(-108.909148353719*nn.^5+ 671.044293885632*nn.^4-1547.0900731642*nn.^3+ ...
1674.14670244375*nn.^2-773.494216929895*nn+130.284211701206)
g1=(coon+i*eee)*e0.*w % 电导率
% Set Ag
% g(20:80,15:J)=g1; % 3.72e8; %0.91*(1e-9)/(36*pi)*w;
%g(1:100,LL:R)=g1; %(-16.4)*(-i)*e0.*w; % 3.72e8; %0.91*(1e-9)/(36*pi)*w;
%e(20:80,80:J)=e0;
% Set Air
for count=1:ArrayNumber
g(ArrayPosX+(count-1)*period:ArrayPosX+(count-1)*period+ArrayWidth,ArrayPosY:ArrayLength)=g1;
end
%k=49;
%for ii=1:1
% g(k:k+2,20:J)=0;
% k=k+5;
%end
%e(20:80,80:J)=e0;
%for k=1:30
% g(i:i+2,140:180)=g1;
% i=i+6;
%e(i:i+1,15:60)=e0;
%end
% -----------------------------------------------------------------------
CA=(1-g*dt/2./e)./(1+g*dt/2./e); % CA=1
CB=dt./e./(1+g*dt/2./e); % CB=0.094248
CP=(1-gm*dt/2/u0)./(1+gm*dt/2/u0); % CP=1
CQ=dt/u0/(1+gm*dt/2/u0); % CQ=6.6315e-007
% 将运算移出循环外，加快程序运行。-----------(Used in Mur boundary condition)
Boundary1=(c*dt-ds)/(c*dt+ds);
Boundary2=c^2*e0*dt/(2*(c*dt+ds));
CORNER=(c*dt-sqrt(2)*ds)./(c*dt+sqrt(2)*ds);
%开始计算
disp('正在计算......');
for n=0:N;
% Switch ----------------------------------------------------------------
if n<pi/w/dt
Switch=0.5*(1-cos(pi*n/(pi/w/dt)));
else
Switch=1;
end
% ------------------------------------------------------------------------
% time evolution ----------------------------------------------------------
x=1:I-1;
y=1:J-1;
Ex1(x,y)=CA(x,y).*Ex(x,y)+CB(x,y).*(Hz(x,y+1)-Hz(x,y))/ds;
Ey1(x,y)=CA(x,y).*Ey(x,y)-CB(x,y).*(Hz(x+1,y)-Hz(x,y))/ds;
x=2:I-1;
y=2:J-1;
Hz1(x,y)=CP*Hz(x,y)-CQ*(Ey1(x,y)-Ey1(x-1,y)-Ex1(x,y)+Ex1(x,y-1))/ds;
Hz1=exp(-i*dt*n*w)*Source*Switch/sqrt(u0/e0)+Hz1;
% -------------------------------------------------------------------------
% 计算截断边界------------------------------------------------------------
%左边界
Hz1(1,y)=Hz(2,y)+Boundary1*(Hz1(2,y)-Hz(1,y))+Boundary2*(Ex1(1,y)-Ex1(1,y-1)+Ex1(2,y)-Ex1(2,y-1));
%右边界
Hz1(I,y)=Hz(I-1,y)-Boundary1*(Hz(I,y)-Hz1(I-1,y))+Boundary2*(Ex1(I-1,y)-Ex1(I-1,y-1)+Ex1(I,y)-Ex1(I,y-1));
%下边界
Hz1(x,1)=Hz(x,2)+Boundary1*(Hz1(x,2)-Hz(x,1))-Boundary2*(Ey1(x,1)-Ey1(x-1,1)+Ey1(x,2)-Ey1(x-1,2));
%上边界
Hz1(x,J)=Hz(x,J-1)-Boundary1*(Hz(x,J)-Hz1(x,J-1))-Boundary2*(Ey1(x,J-1)-Ey1(x-1,J-1)+Ey1(x,J)-Ey1(x-1,J));
%左下角
Hz1(1,1)=Hz(2,2)+CORNER*(Hz1(2,2)-Hz(1,1));
%左上角
Hz1(1,J)=Hz(2,J-1)+CORNER*(Hz1(2,J-1)-Hz(1,J));
%右上角
Hz1(I,J)=Hz(I-1,J-1)-CORNER*(Hz(I,J)-Hz1(I-1,J-1));
%右下角
Hz1(I,1)=Hz(I-1,2)-CORNER*(Hz(I,1)-Hz1(I-1,2));
% -------------------------------------------------------------------------
Hz=Hz1;
Ex=Ex1;
Ey=Ey1;
if mod(n,100)==0
fprintf('the %d step is finished\n',n);
pcolor(0:ds:ds*(J-1),0:ds:ds*(I-1),abs(Hz));
%E=abs(Ex(1:I,1:J)).^2+abs(Ey(1:I,1:J)).^2;
%E=abs(Ex).^2+abs(Ey).^2;
%pcolor(E);
axis equal
shading interp
%shading flat
colorbar('horiz')
xlabel('Y')
ylabel('X')
%pause(0.01)
end
end
%COXIAL
clear
%***********************************************************************
% Fundamental constants
%***********************************************************************
cc=2.99792458e8; %speed of light in free space
muz=4.0*pi*1.0e-7; %permeability of free space
epsz=1.0/(cc*cc*muz); %permittivity of free space
ii=sqrt(-1)
%***********************************************************************
% Grid parameters
%***********************************************************************
ie=34; %number of grid cells in x-direction r=17mm rout=34mm
je=20;
ke=15; %number of grid cells in z-direction l=7.5mm
ib=ie+1;
jb=je+1;
kb=ke+1;
dtheta=2*pi/je;
dx=0.0005; %space increment of cubic lattice 0.5mm
dt=dx/(2.0*cc); %time step
nmax=2000; %total number of time steps
cb=dt/epsz;
db=dt/muz/dx;
er=zeros(ie,jb,kb);
et=zeros(ib,jb,kb);
ez=zeros(ib,jb,ke);
hr=zeros(ib,jb,ke);
ht=zeros(ie,jb,ke);
hz=zeros(ie,jb,kb);
r=zeros(ib,jb,kb);
ez0=zeros(ib,jb);
Nf=1000
fmax=100e+9
df=fmax/Nf
f0=30e+9
ndecay=1/(2*fmax*dt);
n0=1/(2*f0*dt);
t=3.123/dx/ie/2;
%compute the r(i)
for i=1:ib
r(i,:,:)=dx*ie+dx*(i-1);
end
for i=1:ib
ez0(i,:)=besselj(0,t*(ie+i-1)*dx)*bessely(0,t*dx*2*ie)-bessely(0,t*(ie+i-1)*dx)*besselj(0,t*2*dx*ie);
end
%***********************************************************************
% BEGIN TIME-STEPPING LOOP
%***********************************************************************
for n=1:nmax
%***********************************************************************
% Update electric fields
%***********************************************************************
er(1:ie,2:jb,2:ke)=er(1:ie,2:jb,2:ke)+...
dt/dtheta/epsz/(r(1:ie,2:jb,2:ke)+dx/2).*(hz(1:ie,2:jb,2:ke)-hz(1:ie,1:je,2:ke))-...
dt/epsz/dx*(ht(1:ie,2:jb,2:ke)-ht(1:ie,2:jb,1:ke-1));
ez(2:ie,2:jb,1:ke)=ez(2:ie,2:jb,1:ke)+...
dt/epsz/dx*(ht(2:ie,2:jb,1:ke)-ht(1:ie-1,2:jb,1:ke))+dt/epsz/2*(ht(2:ie,2:jb,1:ke)+ht(1:ie-1,2:jb,1:ke))./r(2:ie,2:jb,1:ke)-...
dt/dtheta/epsz*(hr(2:ie,2:jb,1:ke)-hr(2:ie,1:je,1:ke))./r(2:ie,2:jb,1:ke);
et(2:ie,1:jb,2:ke)=et(2:ie,1:jb,2:ke)+dt/epsz/dx*(hr(2:ie,1:jb,2:ke)-hr(2:ie,1:jb,1:ke-1)-hz(2:ie,1:jb,2:ke)+hz(1:ie-1,1:jb,2:ke));
%add excite**************************depend upon the mode which compute.
ez(5,5,3)=exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(5,5,3);
%ez(:,:,3)=ez0(:,:)*exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(:,:,3);
%************* periodic boundary condition for theta direction *****************************
er(1:ie,1,2:ke)=er(1:ie,1,2:ke)+...
dt/dtheta/epsz/(r(1:ie,1,2:ke)+dx/2).*(hz(1:ie,1,2:ke)-hz(1:ie,je,2:ke))-...
dt/epsz/dx*(ht(1:ie,1,2:ke)-ht(1:ie,1,1:ke-1));
ez(2:ie,1,1:ke)=ez(2:ie,1,1:ke)+...
dt/epsz/dx*(ht(2:ie,1,1:ke)-ht(1:ie-1,1,1:ke))+dt/epsz/2*(ht(2:ie,1,1:ke)+ht(1:ie-1,1,1:ke))./r(2:ie,1,1:ke)-...
dt/dtheta/epsz*(hr(2:ie,1,1:ke)-hr(2:ie,je,1:ke))./r(2:ie,1,1:ke);
%***********************************************************************
% Update magnetic fields
%***********************************************************************
hr(2:ie,1:je,1:ke)=hr(2:ie,1:je,1:ke)-...
dt/muz/dtheta*(ez(2:ie,2:jb,1:ke)-ez(2:ie,1:je,1:ke))./r(2:ie,1:je,1:ke)+dt/muz/dx*(et(2:ie,1:je,2:kb)-et(2:ie,1:je,1:ke));
ht(1:ie,1:jb,1:ke)=ht(1:ie,1:jb,1:ke)+...
dt/muz/dx*(er(1:ie,1:jb,1:ke)-er(1:ie,1:jb,2:kb)+ez(2:ib,1:jb,1:ke)-ez(1:ie,1:jb,1:ke));
hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)-dt/muz/2/(r(1:ie,1:je,2:ke)+dx/2).*(et(1:ie,1:je,2:ke)+et(2:ib,1:je,2:ke))-...
dt/muz/dx*(et(2:ib,1:je,2:ke)-et(1:ie,1:je,2:ke))+dt/muz/dtheta/(r(1:ie,1:je,2:ke)+dx/2).*(er(1:ie,2:jb,2:ke)-er(1:ie,1:je,2:ke));
%************* periodic boundary condition for theta direction *****************************
hr(2:ie,jb,1:ke)=hr(2:ie,jb,1:ke)-...
dt/muz/dtheta*(ez(2:ie,2,1:ke)-ez(2:ie,jb,1:ke))./r(2:ie,jb,1:ke)+dt/muz/dx*(et(2:ie,jb,2:kb)-et(2:ie,jb,1:ke));
hz(1:ie,jb,2:ke)=hz(1:ie,jb,2:ke)-dt/muz/2/(r(1:ie,jb,2:ke)+dx/2).*(et(1:ie,jb,2:ke)+et(2:ib,jb,2:ke))-...
dt/muz/dx*(et(2:ib,jb,2:ke)-et(1:ie,jb,2:ke))+dt/muz/dtheta/(r(1:ie,jb,2:ke)+dx/2).*(er(1:ie,2,2:ke)-er(1:ie,jb,2:ke));
%***********************************************************************
%% restroe the imformation at the picked point
%***********************************************************************
gVer1(n)=ez(10,9,10);
end
ff=abs(gVer1);
%plot(ff)
%COXIAL
clear
%***********************************************************************
% Fundamental constants
%***********************************************************************
cc=2.99792458e8; %speed of light in free space
muz=4.0*pi*1.0e-7; %permeability of free space
epsz=1.0/(cc*cc*muz); %permittivity of free space
ii=sqrt(-1)
%***********************************************************************
% Grid parameters
%***********************************************************************
ie=34; %number of grid cells in x-direction r=17mm rout=34mm
je=20;
ke=15; %number of grid cells in z-direction l=7.5mm
ib=ie+1;
jb=je+1;
kb=ke+1;
dtheta=2*pi/je;
dx=0.0005; %space increment of cubic lattice 0.5mm
dt=dx/(2.0*cc); %time step
nmax=2000; %total number of time steps
cb=dt/epsz;
db=dt/muz/dx;
er=zeros(ie,jb,kb);
et=zeros(ib,jb,kb);
ez=zeros(ib,jb,ke);
hr=zeros(ib,jb,ke);
ht=zeros(ie,jb,ke);
hz=zeros(ie,jb,kb);
r=zeros(ib,jb,kb);
ez0=zeros(ib,jb);
Nf=1000
fmax=100e+9
df=fmax/Nf
f0=30e+9
ndecay=1/(2*fmax*dt);
n0=1/(2*f0*dt);
t=3.123/dx/ie/2;
%compute the r(i)
for i=1:ib
r(i,:,:)=dx*ie+dx*(i-1);
end
for i=1:ib
ez0(i,:)=besselj(0,t*(ie+i-1)*dx)*bessely(0,t*dx*2*ie)-bessely(0,t*(ie+i-1)*dx)*besselj(0,t*2*dx*ie);
end
%***********************************************************************
% BEGIN TIME-STEPPING LOOP
%***********************************************************************
for n=1:nmax
%***********************************************************************
% Update electric fields
%***********************************************************************
er(1:ie,2:jb,2:ke)=er(1:ie,2:jb,2:ke)+...
dt/dtheta/epsz/(r(1:ie,2:jb,2:ke)+dx/2).*(hz(1:ie,2:jb,2:ke)-hz(1:ie,1:je,2:ke))-...
dt/epsz/dx*(ht(1:ie,2:jb,2:ke)-ht(1:ie,2:jb,1:ke-1));
ez(2:ie,2:jb,1:ke)=ez(2:ie,2:jb,1:ke)+...
dt/epsz/dx*(ht(2:ie,2:jb,1:ke)-ht(1:ie-1,2:jb,1:ke))+dt/epsz/2*(ht(2:ie,2:jb,1:ke)+ht(1:ie-1,2:jb,1:ke))./r(2:ie,2:jb,1:ke)-...
dt/dtheta/epsz*(hr(2:ie,2:jb,1:ke)-hr(2:ie,1:je,1:ke))./r(2:ie,2:jb,1:ke);
et(2:ie,1:jb,2:ke)=et(2:ie,1:jb,2:ke)+dt/epsz/dx*(hr(2:ie,1:jb,2:ke)-hr(2:ie,1:jb,1:ke-1)-hz(2:ie,1:jb,2:ke)+hz(1:ie-1,1:jb,2:ke));
%add excite**************************depend upon the mode which compute.
ez(5,5,3)=exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(5,5,3);
%ez(:,:,3)=ez0(:,:)*exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(:,:,3);
%************* periodic boundary condition for theta direction *****************************
er(1:ie,1,2:ke)=er(1:ie,1,2:ke)+...
dt/dtheta/epsz/(r(1:ie,1,2:ke)+dx/2).*(hz(1:ie,1,2:ke)-hz(1:ie,je,2:ke))-...
dt/epsz/dx*(ht(1:ie,1,2:ke)-ht(1:ie,1,1:ke-1));
ez(2:ie,1,1:ke)=ez(2:ie,1,1:ke)+...
dt/epsz/dx*(ht(2:ie,1,1:ke)-ht(1:ie-1,1,1:ke))+dt/epsz/2*(ht(2:ie,1,1:ke)+ht(1:ie-1,1,1:ke))./r(2:ie,1,1:ke)-...
dt/dtheta/epsz*(hr(2:ie,1,1:ke)-hr(2:ie,je,1:ke))./r(2:ie,1,1:ke);
%***********************************************************************
% Update magnetic fields
%***********************************************************************
hr(2:ie,1:je,1:ke)=hr(2:ie,1:je,1:ke)-...
dt/muz/dtheta*(ez(2:ie,2:jb,1:ke)-ez(2:ie,1:je,1:ke))./r(2:ie,1:je,1:ke)+dt/muz/dx*(et(2:ie,1:je,2:kb)-et(2:ie,1:je,1:ke));
ht(1:ie,1:jb,1:ke)=ht(1:ie,1:jb,1:ke)+...
dt/muz/dx*(er(1:ie,1:jb,1:ke)-er(1:ie,1:jb,2:kb)+ez(2:ib,1:jb,1:ke)-ez(1:ie,1:jb,1:ke));
hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)-dt/muz/2/(r(1:ie,1:je,2:ke)+dx/2).*(et(1:ie,1:je,2:ke)+et(2:ib,1:je,2:ke))-...
dt/muz/dx*(et(2:ib,1:je,2:ke)-et(1:ie,1:je,2:ke))+dt/muz/dtheta/(r(1:ie,1:je,2:ke)+dx/2).*(er(1:ie,2:jb,2:ke)-er(1:ie,1:je,2:ke));
%************* periodic boundary condition for theta direction *****************************
hr(2:ie,jb,1:ke)=hr(2:ie,jb,1:ke)-...
dt/muz/dtheta*(ez(2:ie,2,1:ke)-ez(2:ie,jb,1:ke))./r(2:ie,jb,1:ke)+dt/muz/dx*(et(2:ie,jb,2:kb)-et(2:ie,jb,1:ke));
hz(1:ie,jb,2:ke)=hz(1:ie,jb,2:ke)-dt/muz/2/(r(1:ie,jb,2:ke)+dx/2).*(et(1:ie,jb,2:ke)+et(2:ib,jb,2:ke))-...
dt/muz/dx*(et(2:ib,jb,2:ke)-et(1:ie,jb,2:ke))+dt/muz/dtheta/(r(1:ie,jb,2:ke)+dx/2).*(er(1:ie,2,2:ke)-er(1:ie,jb,2:ke));
%***********************************************************************
%% restroe the imformation at the picked point
%***********************************************************************
gVer1(n)=ez(10,9,10);
end
ff=abs(gVer1);
%plot(ff)
我把楼主的程序贴出来，大家需要下载的可以到一楼下载
辛苦了，不错
给一部分懒人提供方便，^_^
如果把数组改成动态数组最好。
唉，啥时候我也能达到编程很牛的地步啊，我编程最差了！
编程也没什么，多磨几次，编程就提高了
希望能如楼上所说多练就可以掌握门道。
：23de 我要编程 呵呵 支持
上面的代码作为参考很不错，但不能完全用他的代码作为自己的研究，这样更累
发现有不少学生原封不动地抄袭上面的代码作为作业提交，这样对自己的学习很不负责任，当以后真正做研究要用这部分的知识时，会后悔当初没有真正来学
谢谢分享！！！！！！！！！！！！！！！！！！！！
这是Dennis Sullivan 的书中的程序吧？
谢谢分享！！！！！！！！！！！！！！！！！！！！！
thank you very much!
前面的C语言好像是Dennis Sullivan 的书中的程序啊,
最后两个差不多的是什么程序?是同轴馈线的吗?
下来看看，正学程序！谢谢楼主
kankan!!!!!!!!!!!!1
感谢楼主分享，也感谢二楼的兄弟！
自己学会就好了
哈哈，现在导师要我们写这些东西，看看很有用
Thanks for sharing!
C语言.rar (6.36 KB) 
