projekt_7.m

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

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

function projekt_7()

clear all;
close all;


%%%%%%%%%%%%%%%% frei w鋒lbare Parameter %%%%%%%%%%%%%%%%%%
% Puls Sinus
anzahl_pulse = 3;
pp_wavelength = 40;
freq = 150e6;  % Frequenz des Sinusgenerators

% Anzahl Zellen
KE = 400;

% Dielektrikum
begin_dielek = 320;
width_dielek = 50;
epsilon_r = 2;
epsilon_1 = 2;
delta_0 = 0.25;
omega_0 = 2*pi*100e6;
alpha = delta_0*omega_0;
beta = omega_0*sqrt(1-(delta_0)^2);
gamma = omega_0 / sqrt(1-(delta_0)^2);

% Position Source
source_pos = 50;
T = 0;
% Sensor
sensor_pos = 80;
sensor_buffer = 30;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%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: Sx 
% 5. Zeile: Sxm1
% 6. Zeile: Sxm2
Field = zeros(6,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
spread = anzahl_pulse/2*pp_wavelength;
t0 = (anzahl_pulse*pp_wavelength);
freq = freq*2;
lamda = c0/freq;
dx = lamda / (pp_wavelength); % 6mm

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

factor_1 = 1 / (epsilon_r*epsilon_0);
factor_2 = 2*exp(-alpha*dt)*cos(beta*dt);
factor_3 = -exp(-2*alpha*dt);
factor_4 = exp(-alpha*dt)*sin(beta*dt)*gamma*dt*epsilon_1*epsilon_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: 50 Zellen (3cm bei 500MHz)
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))];
dielek_4 = [zeros(1,begin_dielek) (ones(1,width_dielek).*factor_4) zeros(1,KE - (begin_dielek + width_dielek))];


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

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));
        Field(4,k) = Field(5,k)*dielek_2(k)+Field(6,k)*dielek_3(k)+Field(3,k)*dielek_4(k);
        Field(6,k) = Field(5,k);
        Field(5,k) = Field(4,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(E-Feld, H-Feld, Sensor, Dielektrikum)
        plot(x_axis, Field(1,:), 'b', x_axis, sensorplot, 'r', x_axis, dielektrikumplot, 'm')
        hold off;
        axis([0 KE -2 2]);      
        pause(0.0001);
    end
   
    %Sensor hinlaufende Welle
    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 r點klaufende Welle
    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