#!/usr/bin/octave -qf
%
% 2014
%


%
% calcE
%
% t time
% r radius (might be an array)
% E_theta E field
%

function E_theta = calcE(r, t)

  % r 

  global parameters;
  
  omega = parameters.omega;
  epsilon = parameters.epsilon;
  lambda = parameters.lambda;
  c = parameters.c;

  roh = 2 * pi() * r / lambda;

  E_theta_1 = ( (1 ./ roh) - (1 ./ roh.^3) ) .* cos(roh .- omega*t);
  E_theta_2 = (1 ./roh.^2) .* sin(roh .- omega*t) * -1.;
  E_koeff = omega^3 / (4 * pi() * epsilon * c^3);
  E_theta = -1 * E_koeff * (E_theta_1+E_theta_2);
  
endfunction

%
% calcH
%
% t time
% r radius (might be an array)
% H_theta H field
%

function H_phi = calcH(r, t)

  % r 

  global parameters;
  
  omega = parameters.omega;
  epsilon = parameters.epsilon;
  lambda = parameters.lambda;
  c = parameters.c;

  roh = 2 * pi() * r / lambda;

  H_phi_1 = (1 ./ roh) .* cos(roh .- omega*t);
  H_phi_2 = (1 ./roh.^2) .* sin(roh .- omega*t) * -1.;
  H_koeff = omega^3 / (4 * pi() * epsilon * c^2);
  H_phi = -1 * H_koeff * (H_phi_1 + H_phi_2);
  
endfunction





%
%
% M A I N
%
%

begin_cputime = time();

global parameters;


parameters.framerate = 20;
parameters.duration = 40;
parameters.cycles = 4;
%parameters.framerate = 20;
%parameters.duration = 1;
%parameters.cycles = 1;
parameters.ylim = 1;
parameters.filename = "animation_hertzscher_dipol";
%parameters.Dimensions = "-S800,450";
parameters.Dimensions = "-S1600,900";
parameters.steps = 50;

parameters.linewidth = 2;
%parameters.plot_E = "r;E-Feld;";
%parameters.plot_H = "b;H-Feld;";
parameters.plot_E = "r;E field;";
parameters.plot_H = "b;H field;";
parameters.plot_decay = "k;1/r;";
%parameters.plot_Radius = "Radius normiert zur Wellenlaenge";
%parameters.plot_Amp = "Amplitude normiert bei Radius einer Wellenlaenge";
parameters.plot_Radius = "Radius normalized to wave length";
parameters.plot_Amp = "Amplitude normalized at radius of one wave length";


parameters.f = 0.1; % Hz
parameters.omega = 2 * pi() * parameters.f;
parameters.epsilon = 8.85418781762e-12; % 
parameters.c = 3e8; % m/s
parameters.lambda = parameters.c / parameters.f;

parameters.width = parameters.lambda * 3.; % 3 * wave length

parameters




r = 0:parameters.lambda/parameters.steps:parameters.width;
r_Norm = r ./ parameters.lambda;

t_max = parameters.cycles/parameters.f
%t_schritt = parameters.f/parameters.framerate;
t_schritt = t_max/(parameters.duration*parameters.framerate);
Counter = 1;

Frames = t_max / t_schritt

E_Norm = parameters.omega^3 / (4 * pi() * parameters.epsilon * parameters.c^3);
E_Norm = E_Norm / (2 * pi());

H_Norm = parameters.omega^3 / (4 * pi() * parameters.epsilon * parameters.c^2);
H_Norm = H_Norm / (2*pi());

for t = 0:t_schritt:t_max
  close();

  E = calcE(r, t);
  E = E / E_Norm;

  H = calcH(r, t);
  H = H / H_Norm;
  
  decay = 1 ./ r_Norm;
  
  plot(r_Norm(1,3:end), E(1,3:end), parameters.plot_E, "linewidth", parameters.linewidth);
  hold("on");
  plot(r_Norm(1,3:end), H(1,3:end), parameters.plot_H, "linewidth", parameters.linewidth);
  plot(r_Norm(1,3:end), decay(1,3:end), parameters.plot_decay, "linewidth", parameters.linewidth);
  plot(r_Norm(1,3:end), -1*decay(1,3:end), "color", "k", "linewidth", parameters.linewidth);
  ylim([-15, 15]);
  xlabel(parameters.plot_Radius);
  ylabel(parameters.plot_Amp);
  set (gca(), "linewidth", parameters.linewidth);

  filename = ["./png/", parameters.filename, "-", num2str(Counter, "%04i"), ".png"]
  print(filename, "-dpng", parameters.Dimensions)
  
  Counter += 1;
endfor

printf("Total cpu time: %f seconds\n", time()-begin_cputime);
printf("Finished\n");