% Zeichnen von Eigenvektoren einer FEM-Biegeeigenschwingungsberechnung
%
%  Input:     xyg   -  Knotenkoordinaten, entweder Matrix mit zwei Spalten (x, y)
%                      oder Matrix mit einer Spalte (x)
%             km    -  Koinzidenzmatrix
%             kr    -  Matrix der Randbedingungen (2 oder 3 Spalten)
%             x     -  Matrix, die spaltenweise die Eigenvektoren enthaelt, entweder
%                      pro Knoten 3 Werte (u, v, phi) oder 2 Werte (v, phi)
%             vec1  -  (optional) vec1-te Eigenvektor (vec1-te Spalte aus x) wird
%                      dargestellt (Default: 1)
%             vec2  -  (optional) vec1-te bis vec2-te Eigenvektor werden
%                      dargestellt (Default: vec2 = vec1)
%             magn  -  Verformungen werden in draweigvecs automatisch so
%                      vergroessert, dass eine sinnvolle Groessenordnung im
%                      Verhaeltnis zu den Systemabmessungen entsteht, diese
%                      Vergroesserung kann durch magn veraendert werden
%                      (Default: 1)

function draweigvecs (xyg , km , kr , x , vec1 , vec2 , magn)

[n     , nvecs] = size (x)   ; 
[nnods , dim  ] = size (xyg) ; 
[nnods , kf   ] = size (kr)  ; 

if nargin < 7  enlfac = 1 ; else enlfac = magn ; end

if nargin < 5
    j1 = 1 ;
    j2 = nvecs ;
else
    j1 = vec1 ;
    if nargin < 6
        j2 = j1 ;
    else
        j2 = vec2 ;
    end
end

xy = zeros(nnods , 2) ;
xmax = 0 ;
xmin = 0 ;
ymax = 0 ;
ymin = 0 ;
for i=1:nnods
    xy(i,1) = xyg(i,1) ;
    if dim > 1   xy(i,2) = xyg(i,2)   ; end
    if xy(i,1) > xmax  xmax = xy(i,1) ; end
    if xy(i,1) < xmin  xmin = xy(i,1) ; end
    if xy(i,2) > ymax  ymax = xy(i,2) ; end
    if xy(i,2) < ymin  ymin = xy(i,2) ; end
end

Lmid = (xmax - xmin + ymax - ymin) / 2 ; 

nvec = j2 - j1 + 1 ;
ivec = 1 ;
m    = round (n/kf + 0.1) ;
uv   = zeros (m , kf) ;

Lq   = 0 ;
[nelems , j] = size (km) ;
for j=1:nelems
    nd1 = km (j,1) ;
    nd2 = km (j,2) ;
    Lq = Lq + sqrt((xy(nd2,1)-xy(nd1,1))^2 + (xy(nd2,2)-xy(nd1,2))^2) ;  
end
Lq = Lq/nelems ;                    % "Mittlere Elementlaenge"

for j=j1:j2
    uvq = 0 ;
    k   = 1 ;
    kk  = 1 ;
    while k <= n
      for i=1:kf
         uv(kk,i) = x(k,j) ;
         if i < kf  uvq = uvq + abs(x(k,j)) ; 
         else       uvq = uvq + Lq*abs(x(k,j)) ; 
         end
         k = k + 1 ;
      end
      kk = kk + 1 ;
    end

    uvq = uvq/n ;                   % "Mittlere Knotenverschiebung"
    fac = enlfac*Lmid/(6*uvq) ;
    subplot (nvec,1,ivec) , drawbeamuv (xy , km , kr , uv , fac)  , ...
                            title (strcat(num2str(j),'. Eigenschwingungsform'));
    ivec = ivec + 1 ;
end