projekt_4and5and7.m

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

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

function projekt_4_5_7()

clear all;
close all;


%%%%%%%%%%%%%%%% frei w鋒lbare Parameter %%%%%%%%%%%%%%%%%%
% Puls Sinus
anzahl_pulse = 3;
pp_wavelength = 20;
% Frequenz des Sinusgenerators
freq = 130e6;
%freq = 119.6e6; %inphase




% Anzahl Zellen
KE = 400;

% Dielektrikum
begin_dielek = 320;
width_dielek = 50;

%-------------------------------------------------
epsilon_r = 2;
sigma = 0.01;
chi = 2;
t_0 = 0.001e-6;

%-------------------------------------------------

epsilon_r_3 = 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 = 30;
T = 0;
% Sensor
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];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projekt 5 (Debye Medium)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;


%-------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projekt 4 (lossy Dielektrikum)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Zeile: D-Feld
% 2. Zeile: H-Feld
% 3. Zeile: E-Hilfsfeld
% 4. Zeile: I 
Field2 = zeros(4,KE);

% F黵 ABC
E_end1_2 = 0;
E_end2_2 = 0;
E_anf_2 = 0;

%-------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projekt 7 (Lorentz)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Zeile: D-Feld
% 2. Zeile: H-Feld
% 3. Zeile: E-Hilfsfeld
% 4. Zeile: Sx 
% 5. Zeile: Sxm1
% 6. Zeile: Sxm2
Field3 = zeros(6,KE);

% F黵 ABC
E_end1_3 = 0;
E_end2_3 = 0;
E_anf_3 = 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); 

% 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);
%-------------------------------------------------
factor_1_2 = 1 / (epsilon_r*epsilon_0 + sigma*dt);
factor_2_2 = sigma*dt;
%-------------------------------------------------
factor_1_3 = 1 / (epsilon_r_3*epsilon_0);
factor_2_3 = 2*exp(-alpha*dt)*cos(beta*dt);%*epsilon_0;
factor_3_3 = -exp(-2*alpha*dt);%*epsilon_0;
factor_4_3 = 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);
Fw2 = zeros(1,N);
Fw3 = zeros(1,N);
Fwr = zeros(1,N);
Fwr2 = zeros(1,N);
Fwr3 = 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))];
%-------------------------------------------------
dielek_1_2 = [ones(1,begin_dielek)./epsilon_0 ones(1,width_dielek).*factor_1_2 ones(1,KE - (begin_dielek + width_dielek))./epsilon_0];
dielek_2_2 = [zeros(1,begin_dielek) (ones(1,width_dielek).*factor_2_2) zeros(1,KE - (begin_dielek + width_dielek))];
%-------------------------------------------------
dielek_1_3 = [ones(1,begin_dielek)./epsilon_0 ones(1,width_dielek).*factor_1_3 ones(1,KE - (begin_dielek + width_dielek))./epsilon_0];
dielek_2_3 = [zeros(1,begin_dielek) (ones(1,width_dielek).*factor_2_3) zeros(1,KE - (begin_dielek + width_dielek))];
dielek_3_3 = [zeros(1,begin_dielek) (ones(1,width_dielek).*factor_3_3) zeros(1,KE - (begin_dielek + width_dielek))];
dielek_4_3 = [zeros(1,begin_dielek) (ones(1,width_dielek).*factor_4_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
    
    %-------------------------------------------------
    % D-Feld berechnen
	for k = 2:KE
        Field2(1,k) = Field2(1,k) + ratio_E*(Field2(2,k-1) - Field2(2,k));
        %-------------------------------------------------
        Field3(1,k) = Field3(1,k) + ratio_E*(Field3(2,k-1) - Field3(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; 
        %-------------------------------------------------
        Field2(1,source_pos) =  pulse; 
        %-------------------------------------------------
        Field3(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;
    
    %-------------------------------------------------
    
     for k = 2:KE
        Field2(3,k) = dielek_1_2(k)*( Field2(1,k) - Field2(4,k) );
        Field2(4,k) = Field2(4,k) + dielek_2_2(k)*Field2(3,k);
	end
    
    % Absorbing Boundary Condition
    Field2(3,KE) = E_end2_2;
    E_end2_2 = E_end1_2;
    E_end1_2 = Field2(3,KE-1);
    
    Field2(3,2) = E_anf_2;
    E_anf_2 = Field2(3,1);
    Field2(3,1) = Field2(3,3);
    
	% H-Feld berechnen
	for k = 1:KE-1
        Field2(2,k) = Field2(2,k) + ratio_H*(Field2(3,k) - Field2(3,k+1));
    end;
    
    %-------------------------------------------------
    
    for k = 2:KE
        Field3(3,k) = dielek_1_3(k)*(Field3(1,k) - Field3(4,k));
        Field3(4,k) = Field3(5,k)*dielek_2_3(k)+Field3(6,k)*dielek_3_3(k)+Field3(3,k)*dielek_4_3(k);
        Field3(6,k) = Field3(5,k);
        Field3(5,k) = Field3(4,k);
	end
    
    % Absorbing Boundary Condition
    Field3(3,KE) = E_end2_3;
    E_end2_3 = E_end1_3;
    E_end1_3 = Field3(3,KE-1);
    
    Field3(3,2) = E_anf_3;
    E_anf_3 = Field3(3,1);
    Field3(3,1) = Field3(3,3);
    
    % H-Feld berechnen
	for k = 1:KE-1
        Field3(2,k) = Field3(2,k) + ratio_H*(Field3(3,k) - Field3(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, Field2(1,:), 'r', x_axis, Field3(1,:), 'g' ,x_axis, sensorplot, 'k', x_axis, dielektrikumplot, 'k')
        hold off;
        axis([0 KE -1.5 1.5]);      
        pause(0.00001);
    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);
            Fw2(w) = Fw2(w) +  Field2(1,sensor_pos)*exp(-j*w*sensorcount/N);
            Fw3(w) = Fw3(w) +  Field3(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);
            Fwr2(w) = Fwr2(w) +  Field2(1,sensor_pos)*exp(-j*w*sensorcount2/N);
            Fwr3(w) = Fwr3(w) +  Field3(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(:)),'b')
hold on
plot(abs(Fwr2(:)/Fw2(:)),'r')
plot(abs(Fwr3(:)/Fw3(:)),'g')
title('Reflection');
return