projekt_5.m

这是一个基于FDTD时域有限差分原理的Radar计算程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projekt Numerische Feldberechung   %
%                                    %
%   Marco Angliker, Remo Huber       %
%                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Debye Medium                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function projekt_5()

clear all;
close all;

%%%%%%%%%%%%%%%% frei w鋒lbare Parameter %%%%%%%%%%%%%%%%%%
% Puls Sinus
anzahl_pulse = 3; % Anzahl Sinus
pp_wavelength = 40; %Punke pro Sinus
freq = 1000e6;  % Frequenz des Sinusgenerators

% Anzahl Zellen
KE = 400;

% Dielektrikum
begin_dielek = 320;
width_dielek = 50;
epsilon_r = 2;
sigma = 0.01;
chi = 2;
t_0 = 1/(freq*2);

% Position Source
source_pos = 30;
T = 0;

% Sensor position
sensor_pos = 80;
sensor_buffer = 20;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Plots Dielek & Sensor
sensorplot = zeros(KE, 1);
sensorplot(sensor_pos, 1) = 1;
dielektrikumplot = [zeros(1,begin_dielek) ones(1,width_dielek) zeros(1,KE - (begin_dielek + width_dielek))];

% Zeitachse
x_axis = [1:KE];

% 1. Zeile: D-Feld
% 2. Zeile: H-Feld
% 3. Zeile: E-Hilfsfeld
% 4. Zeile: I 
% 5. Zeile: S
Field = zeros(5,KE);
% F黵 ABC
E_end1 = 0;
E_end2 = 0;
E_anf = 0;

%Konstanten
mue_0 = 1.2566371e-6;
epsilon_0 = 8.85419e-12;
c0 = 299792458;

% Puls Gauss, der 黚er den Sinus 黚erlagert ist
spread = anzahl_pulse/2*pp_wavelength;
t0 = (anzahl_pulse*pp_wavelength);
freq = freq*2; %2 Schritte pro dx
lamda = c0/freq;
dx = lamda / (pp_wavelength); 

% Das 2 im Nennet wird gebraucht, um
% Stabilit鋞sbedingung zu erf黮len.
% 2 Zeitschritte pro Zellenupdate....
dt = dx / (2*c0); 

factor_1 = 1 / (epsilon_r*epsilon_0 + sigma*dt + epsilon_0*chi*(dt/t_0));
factor_2 = sigma*dt;
factor_3 = epsilon_0*chi*(dt/t_0);

ratio_E = dt/(dx);
ratio_H = dt/(mue_0*dx);

figure;
axis manual;
hold on;

% Sensor (Fourieranalyse)

sensorcount = 1;
sensorcount2 = 1;
sensorindex = 1;
N = 2*anzahl_pulse*pp_wavelength;
Fw = zeros(1,N);
Fwr = zeros(1,N);
Fw_temp = zeros(1,N);
Fw_temp2 = zeros(1,N);

% Dielektrikum
dielek_1 = [ones(1,begin_dielek)./epsilon_0 ones(1,width_dielek).*factor_1 ones(1,KE - (begin_dielek + width_dielek))./epsilon_0];
dielek_2 = [zeros(1,begin_dielek) (ones(1,width_dielek).*factor_2) zeros(1,KE - (begin_dielek + width_dielek))];
dielek_3 = [zeros(1,begin_dielek) (ones(1,width_dielek).*factor_3) zeros(1,KE - (begin_dielek + width_dielek))];

% Anzahl Zeitschritte
Nsteps = 2*(begin_dielek-source_pos+begin_dielek+anzahl_pulse*pp_wavelength)

for n = 1:Nsteps

	% D-Feld berechnen
	for k = 2:KE
        Field(1,k) = Field(1,k) + ratio_E*(Field(2,k-1) - Field(2,k));
	end
	
   % *2 in if-Bedingung wegen Stabilit鋞sbdeingung 
   % Siehe -> dt = dx / (2*c0)
   if T <= pp_wavelength*2*anzahl_pulse
        pulse = exp(-1*((T-t0)/spread)^2)*sin(2*pi*freq*dt*T);
        Field(1,source_pos) =  pulse; 
    end;
    
    for k = 2:KE
        Field(3,k) = dielek_1(k)*(Field(1,k) - Field(4,k) - exp(-dt/t_0)*Field(5,k));
        Field(4,k) = Field(4,k) + dielek_2(k)*Field(3,k);
        Field(5,k) = dielek_3(k)*Field(3,k) + exp(-dt/t_0)*Field(5,k);
	end
    
    % Absorbing Boundary Condition
    Field(3,KE) = E_end2;
    E_end2 = E_end1;
    E_end1 = Field(3,KE-1);
    
    Field(3,2) = E_anf;
    E_anf = Field(3,1);
    Field(3,1) = Field(3,3);
    
	% H-Feld berechnen
	for k = 1:KE-1
        Field(2,k) = Field(2,k) + ratio_H*(Field(3,k) - Field(3,k+1));
    end;
    
        % Jede 5. Berechnung der Felder plotten
      if mod(T,5)==0
        clf; % Clear Current Figure
        hold on;
        plot(x_axis, Field(1,:), 'b', x_axis, sensorplot, 'r', x_axis, dielektrikumplot, 'm')
        hold off;
        axis([0 KE -2 2]);      
        pause(0.001);
    end
   
    if (T>=(2*(sensor_pos - source_pos)-sensor_buffer) & T<= 2*((sensor_pos - source_pos) +  anzahl_pulse*pp_wavelength) + sensor_buffer)
        for w = 1:N
            Fw(w) = Fw(w) +  Field(1,sensor_pos)*exp(-j*w*sensorcount/N);
        end
        Fw_temp(sensorcount) = Field(1,sensor_pos);
        sensorcount = sensorcount + 1;
    end
  
    if ( T>= 2*((sensor_pos - source_pos) + 2*(begin_dielek - sensor_pos)) - sensor_buffer & T<= 2*(((sensor_pos - source_pos) + 2*(begin_dielek - sensor_pos)) + anzahl_pulse*pp_wavelength) + sensor_buffer )
        
        for w = 1:N
            Fwr(w) = Fwr(w) +  Field(1,sensor_pos)*exp(-j*w*sensorcount2/N);
        end
        Fw_temp2(sensorcount2) = Field(1,sensor_pos);
        sensorcount2 = sensorcount2 + 1;
        Field(1,sensor_pos);
    end
    
    T = T+1;
end
figure(2)
plot(abs(Fwr(:)/Fw(:)))
title('Reflection');

figure(3);
plot(abs(Fw),'b');
hold on
plot(abs(Fwr),'r');
title('wave spectrum');

figure(4)
plot(Fw_temp,'b');
hold on
plot(Fw_temp2,'r');
title('incident/reflected wave');
return