mu_0   = 1.2566370614*10^-6;
I_0    = 2;
Lambda = 0.69;
omega  = 2*pi*433000000;
t_vals = [0:3.9*10^-11:((pi*4)/2)/omega];
c      = 299792458;

valcount = 401;

x_range = [0:0.05:20];
y_range = zeros(1,valcount);
z_range = zeros(1,valcount);%[-5:0.5:5];

[x_values, y_values, z_values] = meshgrid(x_range, y_range, z_range);
r_values = (x_values.*x_values + y_values.*y_values + z_values.*z_values).^(1/2);

i=0;
for t = t_vals
    i=i+1;
    E_x = zeros(valcount,valcount,valcount);
    E_y = zeros(valcount,valcount,valcount);
    E_z = ((c*c*mu_0*I_0*Lambda)/(4*pi*pi)) ...
          .*((sin((omega/c).*r_values)-sin((omega/c).*r_values-omega*t))./(c.*r_values.^2) ...
          +(cos((omega/c).*r_values-omega*t)-cos((omega/c).*r_values))./(omega.*r_values.^3)) ...
          - ((mu_0*I_0*Lambda*omega)./(4*pi*pi.*r_values)) ...
          .* cos(omega*t-(omega/c).*r_values);

    B_x = ((mu_0*I_0*Lambda)/(4*pi*pi)).*y_values.*((sin((omega/c).*r_values-omega*t)./(r_values.^3))-((omega.*cos((omega/c).*r_values-omega*t))./(c.*r_values.^2)));
    B_y = -((mu_0*I_0*Lambda)/(4*pi*pi)).*x_values.*((sin((omega/c).*r_values-omega*t)./(r_values.^3))-((omega.*cos((omega/c).*r_values-omega*t))./(c.*r_values.^2)));
    B_z = zeros(valcount,valcount,valcount);

    Poynting_x = (-1/mu_0).*E_z.*B_y;
    Poynting_y = (1/mu_0).*E_z.*B_x;
    Poynting_z = zeros(valcount,valcount,valcount);
    %quiver(x_values, y_values, (-1/mu_0).*E_z.*B_y, (1/mu_0).*E_z.*B_x,40);
    %hold on;
    %quiver3(x_values, y_values, z_values, B_x, B_y, B_z);
    %quiver3(x_values, y_values, z_values, E_x, E_y, E_z);
    Poynting_mag = sqrt((Poynting_x.*Poynting_x)+(Poynting_y.*Poynting_y));
    B_mag = sqrt((B_x.*B_x)+(B_y.*B_y));
    plot(x_values(1,:,1), E_z(1,:,1));
    plot(x_values(1,:,1), B_mag(1,:,1));
    xlim([0 20]);
    ylim([-100 100]);
    print(strcat('/home/bram/Desktop/output/', num2str(i)),'-dpng');
end

%pbaspect([1,1,1]);
%title('Poynting Vectors magnitude');
%xlim([-3 3]);
%ylim([-3 3]);
%xlabel('r [m]');
%ylabel('Magnitude poynting vector [W/m^2]');