針對二階精度的時域有限差分程序.

　　現(xiàn)可直接調(diào)用的源信號是:一個周期的正弦信號,高期脈沖,ricker子波.

　　其它信號可手動修改源信號接口,或源生成函數(shù).

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

　　請函數(shù).

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

　　% 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: [email protected]

　　% 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('----------------------------------------------------------------- ');

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

　　子函數(shù)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 step for numerical stability

　　dtmax = 6/7*sqrt(epmin*mumin/(1/dx^2 + 1/dz^2));

　　子函數(shù)2

　　function [dxmax,wlmin,fmax] = finddx(epmax,mumax,srcpulse,t,thres);

　　% finddx.m

　　%

　　% This function finds the maximum spatial discretization that can be used in the

　　% 2-D FDTD modeling codes TM_model2d.m and TE_model2d.m, such that numerical

　　% dispersion is avoided. Second-order accurate time and fourth-order-accurate

　　% spatial derivatives are assumed (i.e., O(2,4)). Consequently, 5 field points

　　% per minimum wavelength are required.

　　%

　　% Note: The dx value obtained with this program is needed to compute the maximum

　　% time step (dt) that can be used to avoid numerical instability. However, the

　　% time vector and source pulse are required in this code to determine the highest

　　% frequency component in the source pulse. For this program, make sure to use a fine