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

　　現可直接調用的源信號是:一個周期的正弦信號,高期脈沖,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 step for numerical stability

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

　　子函數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

　　% temporal discretization for the source pulse, such that no frequency components

　　% present in the pulse are aliased.

　　%

　　% Syntax: [dx,wlmin,fmax] = finddx(epmax,mumax,srcpulse,t,thres)

　　%

　　% where dxmax = maximum spatial discretization possible (m)

　　% wlmin = minimum wavelength in the model (m)

　　% fmax = maximum frequency contained in source pulse (Hz)

　　% epmax = maximum relative dielectric permittivity in grid

　　% mumax = maximum relative magnetic permeability in grid

　　% srcpulse = source pulse for FDTD simulation

　　% t = associated time vector (s)

　　% thres = threshold to determine maximum frequency in source pulse

　　% (default = 0.02)

　　%

　　% by James Irving

　　% July 2005

　　if nargin==4; thres=0.02; end

　　% convert relative permittivity and permeability to true values

　　mu0 = 1.2566370614e-6;

　　ep0 = 8.8541878176e-12;

　　epmax = epmax*ep0;

　　mumax = mumax*mu0;

　　% compute amplitude spectrum of source pulse and corresponding frequency vector

　　n = 2^nextpow2(length(srcpulse));

　　W = abs(fftshift(fft(srcpulse,n)));

　　W = W./max(W);

　　fn = 0.5/(t(2)-t(1));

　　df = 2.*fn/n;

　　f = -fn:df:fn-df;

　　W = W(n/2+1:end);

　　f = f(n/2+1:end);

　　% determine the maximum allowable spatial disretization

　　% (5 grid points per minimum wavelength are needed to avoid dispersion)

　　fmax = f(max(find(W>=thres)));

　　wlmin = 1/(fmax*sqrt(epmax*mumax));

　　dxmax = wlmin/5;

　　子函數3

　　function [excitation]=gprmaxso(type,amp,freq,dt,total_time);

　　% GPRMAXSO Computes the excitation function used in 'GprMax2D/3D'

　　% simulators for ground probing radar.

　　%

　　% [excitation] = gprmaxso('source_type',Amplitude,frequency,Time_step,Time_window)

　　% source_type can be 'cont_sine', 'gaussian', 'ricker'

　　% Amplitude is the amplitude of the source

　　% frequency is the frequency of the source in Hz

　　% Time_step is the time step in seconds

　　% Time_window is the total simulated time in seconds

　　%

　　% excitation is a vector which contains the excitation function.

　　% If you type plot(excitation) Matlab will plot it.

　　% You can use the signal processing capabilities of Matlab

　　% to get a Spectrum etc.

　　%

　　% Copyright: Antonis Giannopoulos, 2002 This file can be distributed freely.

　　RAMPD=0.25;

　　if(nargin < 5)

　　error('GPRMAXSO requires all five arguments ');

　　end;

　　if(isstr(type)~=1)

　　error('First argument should be a source type');

　　end;

　　if(freq==0)

　　error(['Frequency is zero']);

　　end;

　　iter=total_time/dt;

　　time=0;

　　if(strcmp(type,'ricker')==1)

　　rickth=2.0*pi*pi*freq*freq;

　　rickper=1.0/freq;

　　ricksc=sqrt(exp(1.0)/(2.0*rickth));

　　i=1;

　　while(time<=total_time)

　　delay=(time-rickper);

　　temp=exp(-rickth*delay*delay);

　　excitation(i)=ricksc*temp*(-2.0)*rickth*delay;

　　time=time+dt;

　　i=i+1;

　　end;

　　end;

　　if(strcmp(type,'gaussian')==1)

　　rickper=1.0/freq;

　　rickth=2.0*pi*pi*freq*freq;

　　i=1;

　　while(time<=total_time)

　　delay=(time-rickper);

　　excitation(i)=exp((-rickth)*delay*delay);

　　time=time+dt;

　　i=i+1;

　　end;

　　end;

　　if(strcmp(type,'cont_sine')==1)

　　i=1;

　　while(time<=total_time)

　　ramp=time*RAMPD*freq;

　　if(ramp>1.0)

　　ramp=1.0;

　　end;

　　excitation(i)=ramp*sin(2.0*pi*freq*time);

　　time=time+dt;

　　i=i+1;

　　end;

　　end;

　　if(strcmp(type,'sine')==1)

　　i=1;

　　while(time<=total_time)

　　excitation(i)=sin(2.0*pi*freq*time);

　　if(time*freq>1.0)

　　excitation(i)=0.0;

　　end;

　　time=time+dt;

　　i=i+1;

　　end;

　　end;

　　excitation=excitation.*amp;