% Liebe Kolleginnen und Kollegen,
%
% ich zeige Ihnen hier, wie man nach der Berechnung einer
% Temperaturverteilung (auf dem kartesischen Netz) die
% Temperaturgradienten durch Pfeile grafisch darstellen
% kann.
% Diese Darstellung ist ausserordentlich brauchbar fuer die
% physikalische Interpretation der Ergebnisse: 


disp('Berechnung der Temp.verteilung auf dem quadr. Netz:');
% Berechnung des Temp-Files fuer Matlab:
  tfile=zeros(101,101);
  for iy=1:101
    y=-0.5+0.04*(iy-1);
    for ix=1:101
      x=0.04*(ix-1);  
  .
  .
  .
      tfile(ix,iy)=c(ind,1)+c(ind,2)*u0+c(ind,3)*v0;
    end
  end

  pcolor(tfile');
  caxis([0,150]);
  axis('equal');
  title('AK NUMERIK FE stat. Waermeverteilung'); 

  hold on

% Hier beginnt die Berechnung und grafische Darstellung der Temperatur-
% gradienten (= des Waermeflusses) auf der Platte:

for iy=2:5:100
  for ix=2:5:100
    xc=0.04*(ix-1); yc=-0.5+0.04*(iy-1);

% Der "Waermepfeil" wird nur gerechnet, wenn unter den relevanten Temperaturen
% keine Null vorkommt:
    T1=tfile(ix+1,iy); T2=tfile(ix-1,iy); T3=tfile(ix,iy+1); T4=tfile(ix,iy-1);
    vTemp=[T1 T2 T3 T4];
    if(sum(vTemp==0)==0)
      dTdx=(T1-T2)/(2*0.04);
      dTdy=(T3-T4)/(2*0.04);
      f=0.0075;       % dieser (willkuerliche) Faktor hat nur eine grafische,
                      % keine physikalische Bedeutung.

      wflussx=-dTdx*f;   % Beachten Sie das negative Vorzeichen: die Waerme 
      wflussy=-dTdy*f;   % fliesst stets in Richtung der fallenden Temperatur!

      ixw=ix+wflussx/0.04;
      iyw=iy+wflussy/0.04;
      arrz(ix,ixw,iy,iyw); % dieses Matlab-Skript, welches den Pfeil zeichnet, 
                           % habe ich verbrochen (s.u.). Ich bin sicher, 
                           % irgendwo im Matlab-Repertoire gibt's i
                           % was Besseres ...
    end
  end
end






%
%
%  function y=arrz(x1,x2,y1,y2)
%
%  L=sqrt((x2-x1)^2+(y2-y1)^2);
%  e1=[x2-x1 y2-y1]/L;
%  P=0.8*L*e1+[x1 y1];
%
%  e2=[y2-y1,x1-x2]/L;
%  S1=P+e2*L/20;
%  S2=P-e2*L/20;
%  plot([x1 x2],[y1 y2],'w-');
%  hold on
%  plot([x2 S1(1) S2(1) x2],[y2 S1(2) S2(2) y2],'w-');


% H.Sormann  Mai 2011