postproc.m

非连续波导时间域有限差分法仿真（TD-FDM）

% post procession
% time domain -> frequency domain
% visulization

% solve reflection wave for
% magnetic field at incoming port
sample_freq = 1/dt;
f = sample_freq * (0:N/2-1)/N;
df = sample_freq / N;

total = Hz(round(Ny/2),:);
trans = Hz(node-(Ny-round(Ny/2)),:);
ref = total - g;

G = fft(g,N);
pG = G.*conj(G) / N;
Ref = fft(ref,N);
pRef = Ref.*conj(Ref) / N;
Trans = fft(trans,N);
pTrans = Trans.*conj(Trans) / N;

refCoeff = pRef./pG;
transCoeff = pTrans./pG;

startIdx = round(300e6/df)-1;
endIdx = round(900e6/df)+1;
figure;
plot(f(startIdx:endIdx), refCoeff(startIdx:endIdx));
%title('Reflection Coefficient versus Frequency');
xlabel('Frequency(Hz)');
ylabel('Reflection coefficient versus Freq');

figure;
plot(f(startIdx:endIdx), transCoeff(startIdx:endIdx));
xlabel('Frequency(Hz)');
ylabel('Transition coefficient versus Freq');

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

figure;
plot(t,g);
title('Gaussian pulse versus Time');
hold on;

figure;
cnt1 = nmbr(round(Ny/2),1);
plot(t,Hz(cnt1,:),'r');
hold on;
%title('Field value at incoming port');
%figure;
cnt2 = nmbr(round(Ny/2),MAX_X);
plot(t,Hz(cnt2,:),'b');
%title('Field value at outgoing port');
legend('H at incomg port','H at outgoing port');
title('Magneitc field at two ports');

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

disp('Total time for time update...');
tt
% make a movie, fun stuff!
makeMovie;