#!/usr/bin/octave 

%
% Diagramm zur Dipolantenne mit Länge zu Eingangsimpedanz
%
% aus Krak "Antennen und Strahlungsfelder" bzw. http://rediscoveries.blogspot.de/2011/08/half-wavelength-dipole-that-isnt.html
%
% 2014
%

% Abhängigkeit: pkg install -auto -install specfun


%
% reaktanz()
%

function [R_E, X_E] = reaktanz(h, a)
global k_0;
global Z_0;
global C;

  Term_0 = Z_0 / (4*pi());
  Term_1 = - Ci(2*k_0*h) + C + log(2*k_0*h);
  Term_2 = Si(4*k_0*h) - 2 * Si(2*k_0*h);
  Term_3 = Ci(4*k_0*h) - 2 * Ci(2*k_0*h) + C + log(k_0*h);
  Term_4 = 1 ./ (sin(k_0 * h)).^2;

  R_E = Term_0 * Term_4 .* (2 * Term_1 + sin(2 * k_0 * h) .* Term_2 + cos(2*k_0*h) .* Term_3);


  Term_X_1 = Si( 2 * k_0 * h);
  Term_X_2 = Term_2 * -1;
  Term_X_3 = Ci(4*k_0*h) - 2 * Ci(2*k_0*h) + C + log((k_0 * a).^2) - log(k_0*h);

  X_E = Term_0 * Term_4 .* (2 * Term_X_1 + cos(2 * k_0 * h) .* Term_X_2 + sin(2*k_0*h) .* Term_X_3);

endfunction


%
% M  A  I  N
%

% Normierung
lambda = 1.;

% Länge der Dipolantenne
schritt = 0.01;
% schritt = 0.025
l = schritt:schritt:3*lambda;
ratio = 20;
%l = 0.5;
% l = 0.5*lambda;
%l
h = l/2;

% Freiraumimpedanz nach WP
global Z_0 = 376.730313461; %  Ohm
% Z_0 = 1;

% Euler-Mascheroni-Konstante nach WP
global C = 0.5772156649;

% Kreiswellenzahl
global k_0 = 2*pi()/lambda; % TODO

% http://essay.utwente.nl/61459/1/MSc_P_Prins.pdf
% eta = Z_0 / Z_m = Z_E / Z_0 = Z_S
%printf("Master\n");
%Term_1 = Ci(2*k_0*h) + C + log(2*k_0*h);
%Term_2 = Si(4*k_0*h) - 2 * Si(2*k_0*h);
%Term_3 = - Ci(4*k_0*h) + 2 * Ci(2*k_0*h) + C + log(k_0*h);
%Term_4 = 1 ./ (sin(k_0 * h)).^2;


%R_E

% Kark [Zuh53]
printf("Kark\n");

a_1000 = h/1000;
[R_E, X_E1000] = reaktanz(h, a_1000);
Z_E1000 = sqrt(R_E.^2 + X_E1000.^2);

a_10 = h/10;
[R_E, X_E10] = reaktanz(h, a_10);
Z_E10 = sqrt(R_E.^2 + X_E10.^2);

close();

% set (0, 'DefaultFigurePosition', [1 0.5 * get(0, "screensize")(4:-1:3) get(0, "DefaultFigurePosition")(4)]);

figure('Position',[100, 100, 1200,1000]);

plot( l/lambda, R_E, "k;R;", "linewidth", 2);
hold("on");
plot( l/lambda, X_E1000, "--r;X_{1000};", "linewidth", 2);
plot(l/lambda, Z_E1000, "r;|Z_{1000}|;", "linewidth", 2);
plot( l/lambda, X_E10, "--b;X_{10};", "linewidth", 2);
plot(l/lambda, Z_E10, "b;|Z_{10}|;", "linewidth", 2);

ylim([-250, 250]);
xlabel("Ratio l/\\lambda");
ylabel("Impedance in Ohm");

legend ("location", "southeast");

print("dipol_impedance.png", "-dpng", "-S1600,900")

% printf("Weiter mit beliebiger Taste...\n");
% pause();