針對二階精度的時域有限差分程序.
現(xiàn)可直接調用的源信號是:一個周期的正弦信號,高期脈沖,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: 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('----------------------------------------------------------------- ');
%**************************************************
子函數(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[!--empirenews.page--]
% 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;
子函數(shù)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)[!--empirenews.page--]
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;
下載文檔
