FDTD参数选择估计程序

针对二阶精度的时域有限差分程序.

现可直接调用的源信号是:一个周期的正弦信号,高期脉冲,ricker子波.

其它信号可手动修改源信号接口,或源生成函数.

---------------

请函数.

%************************************************************

% 1. determine maximum possible spatial field discretization.

% (in order to avoid numerical dispersion).(5 grid points per

% minimum wavelength are needed to avoid dispersion).

% 2. find the maximum possible time step using this dx and dz.

% (in order to avoid numerical instability).

% Coded by yiling. Email: yiling@email.jlu.edu.cn

% Date: 2008

%*************************************************************************+

clear;

clc;

%--------------------------------------------------------------------------

dx=0.02; % (m)

dy=0.02; % (m)

epsilonmax=25; % Epsion. maximum relative dielectric permittivity.

mumax=1; % Mu. maximum relative magnetic permeability.

sourcetype='ricker'; % can be 'cont_sine', 'gaussian', 'ricker'.

freq=100e6; % (Hz)

amp=1; % amplitude.

thres=0.02; % threshold to determine maximum frequency in source pulse.(proposed = 0.02).

%--------------------------------------------------------------------------

Timewindows=528; % (ns)

%--------------------------------------------------------------------------

%*************************************************************************+

%--------------------------------------------------------------------------

vlight=0.3;

epsilonmin=1; % Epsion. minimum relative dielectric permittivity.

mumin=1; % Mu. minimum relative magnetic permeability.

%--------------------------------------------------------------------------

dt=1/(vlight*sqrt(1/dx^2+1/dy^2));

% minwavelength=vlight/sqrt(epsilinmax);

%--------------------------------------------------------------------------

t=0:dt:Timewindows;

dt=dt*1e-9;

t=t*1e-9;

Timewindows=Timewindows*1e-9;

source=gprmaxso(sourcetype,amp,freq,dt,Timewindows);

[dxmax,wlmin,fmax] = finddx(epsilonmax,mumax,source,t,thres);

%--------------------------------------------------------------------------

disp('----------------------------------------------------------------- ');

disp(['Maximum frequency contained in source pulse = ',num2str(fmax/1e6),' MHz']);

disp(['Minimum wavelength in simulation grid = ',num2str(wlmin),' m']);

disp(['Maximum possible electric/magnetic field discretization (dx,dy) = ',num2str(dxmax),' m']);

disp(' ');

%--------------------------------------------------------------------------

%--------------------------------------------------------------------------

dtmax = finddt(epsilonmin,mumin,dxmax,dxmax);

disp(['Maximum possible time step with this discretization = ',num2str(dtmax/1e-9),' ns']);

disp('----------------------------------------------------------------- ');

%**************************************************

子函数1

function dtmax = finddt(epmin,mumin,dx,dz);

% finddt.m

%

% This function finds the maximum time step that can be used in the 2-D

% FDTD modeling codes TM_model2d.m and TE_model2d.m, such that they remain

% numerically stable. Second-order-accurate time and fourth-order-accurate

% spatial derivatives are assumed (i.e., O(2,4)).

%

% Syntax: dtmax = finddt(epmin,mumin,dx,dz)

%

% where dtmax = maximum time step for FDTD to be stable

% epmin = minimum relative dielectric permittivity in grid

% mumin = minimum relative magnetic permeability in grid

% dx = spatial discretization in x-direction (m)

% dz = spatial discretization in z-direction (m)

%

% by James Irving

% July 2005

% convert relative permittivity and permeability to true values

mu0 = 1.2566370614e-6;

ep0 = 8.8541878176e-12;

epmin = epmin*ep0;

mumin = mumin*mu0;

% determine maximum allowable time s