projekt_3.m

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

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

function projekt_3()

clear all;
close all;

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

% Anzahl Zellen
KE = 400;

% Dielektrikum
begin_dielek = 320;
width_dielek = 50;
epsilon_r = 4;
sigma = 0.04;

% 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: E-Feld
% 2. Zeile: H-Feld
Field = zeros(2,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_E1 = (1-(dt*sigma/(2*epsilon_r*epsilon_0)))/(1+(dt*sigma/(2*epsilon_r*epsilon_0)));
factor_E2 = 1/(1+(dt*sigma/(2*epsilon_r*epsilon_0)));
ratio_E = dt/(epsilon_0*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: 50 Zellen (3cm bei 500MHz)
dielek_E = [ones(1,begin_dielek) (ones(1,width_dielek).*1/epsilon_r).*factor_E2 ones(1,KE - (begin_dielek + width_dielek))];
dielek_E2 = [ones(1,begin_dielek) (ones(1,width_dielek).*factor_E1) ones(1,KE - (begin_dielek + width_dielek))];

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

for n = 1:Nsteps

    % E-Feld berechnen
	for k = 2:KE
        % ex[k] = ex[k] + 0.5*(hy[k-1] - h[k])
        Field(1,k) = dielek_E2(k)*Field(1,k) + dielek_E(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;
    
    % Absorbing Boundary Condition
    Field(1,KE) = E_end2;
    E_end2 = E_end1;
    E_end1 = Field(1,KE-1);
    
    Field(1,2) = E_anf;
    E_anf = Field(1,1); %Field(1,1) wird hier als Buffer gebraucht
    Field(1,1) = Field(1,3);
	
	% H-Feld berechnen
	for k = 1:KE-1
        % hy[k] = hy[k] + 0.5*(ex[k] - ex[k+1]
        Field(2,k) = Field(2,k) + ratio_H*(Field(1,k) - Field(1,k+1));
    end;
        % Jede 5. Berechnung der Felder plotten
      if mod(T,5)==0
        clf; % Clear Current Figure
        hold on;
        %H-Feld auf dem Plot skaliert, damit man es sieht
        plot(x_axis, Field(1,:), 'b', x_axis, 1.4193e+002*Field(2,:), 'g', x_axis, sensorplot, 'r', x_axis, dielektrikumplot, 'm')
        hold off;
        axis([0 KE -2 2]);      
        pause(0.001);
    end
   
    %Sensor first time
    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
  
    %Sensor seconde time
    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

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