#! /usr/bin/octave -q
% ---------------------------------------------------------
% Dieses Program laeuft unter GNU-Octave
% Es berechnet die Streuung einer ebenen Welle an
% einer Kugel
%
% Klaus.Ehrenfried@dlr.de
%

% ---------------------------------------------------------
% Invertiere Diagonalmatrix; wird zur Berechnung der
% Pseudo-Inversen benoetigt
%

function SI = invert_diagonal(S)

TINY=1e-10;
[m,n] = size(S);
SI = zeros([n,m]);

for i=1:min(m,n)
  if (abs(S(i,i)) > TINY)
    SI(i,i) = 1.0/S(i,i);
  end
end

endfunction

% -------------------------------------------------------
% Berechne Staerke und Ableitung der Felder von Einheits-Monopolen 
% in Normalenrichtung
% (tx,ty,tz): Beobachtungs/Testpunkt
% (nx,ny,nz): Normalenvektor
% (qx,qy,qy): Koordinaten der Quellen; qx,qy,qz sind Vektoren;
%             length(qx) ist die Anzahl der Quellen
% k_welle: Wellenzahl
%

function [phi, dphi_dn] = ber_monopol(tx,ty,tz,nx,ny,nz,qx,qy,qz,k_welle)
dx = tx - qx;
dy = ty - qy;
dz = tz - qz;
r = sqrt(dx.*dx + dy.*dy + dz.*dz);
dr_dn = (dx * nx + dy * ny + dz * nz) ./ r;
dphi_dr = (-1 ./ (r .* r) - sqrt(-1.0) * k_welle ./ r) ...
	  .* exp(-sqrt(-1.0) * k_welle * r);
dphi_dn = dphi_dr .* dr_dn;
phi = 1 ./ r .* exp(-sqrt(-1.0)*k_welle*r);
endfunction

% -------------------------------------------------------
% Matrix berechnen
% A_ij = dphi_dn am Punkt (i) von Einheitsquelle am Punkt (j) verursacht
% (tx,ty,tz): Testpunkte
% (nx,ny,nz): Normelenvektoren in Testpunkten
% (qx,qy,qz): Quellpostionen
% Nt: Anzahl der Testpunkte; ist gleich length(tx)
% Nq: Anzahl der Quellen; ist gleich length(qx)
% k_welle: Wellenzahl
%

function A = berechne_matrix(tx,ty,tz,nx,ny,nz,qx,qy,qz,Nt,Nq,k_welle)
A = zeros(Nt,Nq);

for i=1:Nt
  if (mod(i,100) == 0)
    printf('loop: %d\n', i);
  end
  [phi, dphi_dn] = ...
      ber_monopol(tx(i),ty(i),tz(i),nx(i),ny(i),nz(i),qx,qy,qz,k_welle);
  A(i,:) = dphi_dn;
end

%disp(abs(A(1:10,1:10)));

printf('Kondition: %e\n', cond(A));
endfunction

% -------------------------------------------------------
% Berechne Ableitung des Feldes in einer ebene Welle der Staerke 1
% in Normalenrichtung
% (tx,ty,tz): Testpunkte
% (nx,ny,nz): Normelenvektoren in Testpunkten
% (nex,ney,nez): Ausbreitungsrichtung der Welle
% k_welle: Wellenzahl
%

function [phi, dphi_dn] = ebene_welle(tx,ty,tz,nx,ny,nz,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));
dphi_dx = -sqrt(-1) * k_x * phi;
dphi_dy = -sqrt(-1) * k_y * phi;
dphi_dz = -sqrt(-1) * k_z * phi;

dphi_dn = (dphi_dx .* nx + dphi_dy .* ny + dphi_dz .* nz);
endfunction

% -------------------------------------------------------
% Berechnete ueberlagerte Loesung aller Quellen und
% deren Normalableitung
% (xp,yp,zp): Testpunkte
% (nxp,nyp,nzp): Normalenvektoren
% k_welle: Wellenzahl
% (qx,qy,qz): Quellpostionen
% q: komplexe Quellstaerke der Quellen
% Ndim: Anzahl der Testpunkte
%

function [pq,vnqp] = ...
    teste_loesung(xp,yp,zp,nxp,nyp,nzp,k_welle,qx,qy,qz,q,Ndim)

vnqp = zeros(1,Ndim);
pq = zeros(1,Ndim);

for i=1:Ndim
  [phiq, dphiq_dn] = ...
      ber_monopol(xp(i),yp(i),zp(i),nxp(i),nyp(i),nzp(i),qx,qy,qz,k_welle);

  vnqp(i) = dphiq_dn * q;
  pq(i) = phiq * q;
end

endfunction

% -------------------------------------------------------
% Erzeugen der Koordinaten fuer Testpunkte
%

Nt = 201;
theta = linspace(0.0,pi,Nt);

tx = cos(theta);
ty = sin(theta);
tz = zeros(1,Nt);

nx = tx;
ny = ty;
nz = tz;

nnr = sqrt(nx.*nx+ny.*ny+nz.*nz);
nx = nx ./ nnr;
ny = ny ./ nnr;
nz = nz ./ nnr;

% -------------------------------------------------------
% Erzeugen der Koordinaten fuer Quellen
%

Nq = 201;

stretch = 0.8;
qx = linspace(-stretch, stretch, Nq);
qy = zeros(1,Nq);
qz = zeros(1,Nq);

xlabel('x');
ylabel('y');
plot(qx,qy,'o3;Quellen;',tx,ty,'-+1;Testpunkte;');
printf('Enter druecken!\n');
pause;

if (length(tx) < length(qx))
  printf('Error: Daten passen nicht zusammen!\n');
  exit;
end

% -------------------------------------------------------
% Ebene Welle
%

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

printf('Ebene Welle: (%f,%f,%f)\n', nex, ney, nez);

% -------------------------------------------------------
% Punkte zur Darstellung der Loseung auf der Oberflaeche
%

Ndim = 501;

thetap = linspace(0.0,pi,Ndim);

xp = cos(thetap);
yp = zeros(1,Ndim);
zp = sin(thetap);

nxp = xp;
nyp = yp;
nzp = zp;

nnrp = sqrt(nxp.*nxp+nyp.*nyp+nzp.*nzp);
nxp = nxp ./ nnrp;
nyp = nyp ./ nnrp;
nzp = nzp ./ nnrp;

% -------------------------------------------------------
% Aufstellen der Matrix
%

%Nt = length(tx);
%Nq = length(qx);

printf('Anzahl der Testpunkte: %d\n', Nt);
printf('Anzahl der Quellpunkte: %d\n', Nq);

fr = 1.7;

lambda = 1.0/fr;
printf('lambda: %f\n', lambda);
k_welle = 2.0 * pi / lambda;
printf('k: %e\n', k_welle);

A = berechne_matrix(tx,ty,tz,nx,ny,nz,qx,qy,qz,Nt,Nq,k_welle);
[pe, vne] = ebene_welle(tx,ty,tz,nx,ny,nz,nex,ney,nez,k_welle);

% -------------------------------------------------------
% Loesung berechnen ueber Pseudo-Inverse von A.
% Falls Nt = Nq ist
% Kondition von A oft relativ hoch, daher Weg ueber
% Pseudo-Inverse auch bei Nt = Nq robuster
%

[U,S,V] = svd(A);
SI = invert_diagonal(S);

q = -V * SI * U' * vne.';

rtest = A * q + vne.';

printf('Maximaler Rest: %e\n', max(abs(rtest)));
printf('Relativer Fehler: %e\n', max(abs(rtest))/max(abs(vne)));
printf('norm. RMS Fehler: %e\n', ...
       sqrt(sum(abs(rtest).**2)/length(rtest))/max(abs(vne)));

[pq, vnqp] = teste_loesung(xp,yp,zp,nxp,nyp,nzp,k_welle,qx,qy,qz,q,Ndim);
[pe, vnep] = ebene_welle(xp,yp,zp,nxp,nyp,nzp,nex,ney,nez,k_welle);

xlabel('x');
ylabel('dphi_dn');
plot(xp,real(vnep),'-o3;Re vne;',xp,-real(vnqp),'-4;-Re vnq;');
printf('Enter druecken!\n');
pause;
plot(xp,imag(vnep),'-o3;Im vne;',xp,-imag(vnqp),'-4;-Im vnq;');
printf('Enter druecken!\n');
pause;
ylabel('phi');
plot(xp,real(pe),'-o2;Re pe;',xp,real(pq),'-6;Re pq;');
printf('Enter druecken!\n');
pause;
ylabel('q');
plot(qx,real(q),'-o2;Re q;',qx,imag(q),'-6;Im q;');
printf('Enter druecken!\n');
pause;

% -------------------------------------------------------
% Ergebnis abspeichern
%

save -ascii loesung_k4.data qx qy qz q k_welle;

unit=fopen('testpunkte.dat','w');
data = [tx;ty;tz];
fprintf(unit,'%e %e %e\n',data);
fclose(unit);

unit=fopen('quellpunkte.dat','w');
data = [qx;qy;qz];
fprintf(unit,'%e %e %e\n',data);
fclose(unit);

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