FDTD參數選擇估計程序

　　% 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)