#! /usr/bin/octave -q
% -------------------------------------------------------
% Dieses Program laeuft unter GNU-Octave
% Es berechnet die beiden Teilloesungen und die Gesamtloesung
% in einer Ebene aus den zuvor berecheten Quellstaerken
%
% Klaus.Ehrenfried@dlr.de
%
%

load loesung_k4.data;
load testpunkte.data;

whos;

% -------------------------------------------------------
% Berechne Felder von Einheits-Monopolen 
% (tx,ty,tz): Beobachtungs/Testpunkt
% (qx,qy,qy): Koordinaten der Quellen; qx,qy,qz sind Vektoren;
%             length(qx) ist die Anzahl der Quellen
% k_welle: Wellenzahl
%

function phi = phi_monopol(tx,ty,tz,qx,qy,qz,k_welle)
dx = tx - qx;
dy = ty - qy;
dz = tz - qz;
r = sqrt(dx.*dx + dy.*dy + dz.*dz);
phi = 1 ./ r .* exp(-sqrt(-1.0)*k_welle*r);
endfunction

% -------------------------------------------------------
% Berechne Feld einer ebene Welle
% (tx,ty,tz): Testpunkte
% (nex,ney,nez): Ausbreitungsrichtung der Welle
% k_welle: Wellenzahl
%

function phi = phi_ebene_welle(tx,ty,tz,nex,ney,nez,k_welle)
k_x = k_welle * nex;
k_y = k_welle * ney;
k_z = k_welle * nez;

phi = exp(-sqrt(-1) * (k_x * tx + k_y * ty + k_z * tz));
endfunction

% -------------------------------------------------------------------
% Richtung der ebenen Welle
%

nex = 1.0;
ney = 0.0;
nez = 0.0;

% -------------------------------------------------------------------
% Gitter erzeugen
%

Nwinkel = 361;
Nradius = 161;
Rmax = 8.0;

walpha = linspace(0.0,2.0*pi,Nwinkel);
rr = linspace(1.0, Rmax, Nradius);

walpha_m = repmat(walpha, Nradius, 1);
rr_m = repmat(rr',1,Nwinkel);

xx = rr_m .* cos(walpha_m);
yy = rr_m .* sin(walpha_m);
zz = zeros(Nradius,Nwinkel);
phi_r = zeros(Nradius,Nwinkel);

printf('Berechne Feld mit %d x %d Punkten...\n', Nradius, Nwinkel);
for i=1:Nradius
  printf('%d/%d\n',i,Nradius);
  for j=1:Nwinkel
    phi_1 = phi_monopol(xx(i,j),yy(i,j),zz(i,j),qx,qy,qz,k_welle);
    phi_r(i,j) = phi_1 * q;
  end
end

phi_e = phi_ebene_welle(xx,yy,zz,nex,ney,nez,k_welle);
phi_ges = phi_e + phi_r;

printf('Abspeichern...\n');

%contour(xx,yy,zz);
%pause;

unit=fopen('feld.data','w');

for n=1:Nradius
  data = ...
      [xx(n,:);yy(n,:);real(phi_e(n,:));real(phi_r(n,:));real(phi_ges(n,:))];
  fprintf(unit,'%e %e %e %e %e\n',data);
  fprintf(unit,'\n');
end

fclose(unit);

% -------------------------------------------------------------------
% Ende
% -------------------------------------------------------------------
