Izvorno besedilo

Računalnik je prislovično hiter, vesten in neumoren delavec, le bister ni kaj dosti. Stori natanko samo tisto, kar smo mu v programu naročili. Da se izognemo naročanju vedno enakih korakov, imajo programski jeziki ukaze, ki omogočajo ponavljanja – zanke. Skoraj v vsaki nalogi doslej smo uporabili tako zanko (<tt>for i := 1 to n do</tt>).

Gibkost programiranja pa zagotavljajo predvsem razvejitveni ukazi, ko je nadaljnji potek računanja odvisen od izpolnitve kontrolnega pogoja (<tt>if '''pogoj''' then do ... else do'' ...''</tt>). Posebna vrsta zelo gibkih razvejitev so klici podprogramov in funkcij. S temi orodji bomo skicirali zgradbo programa za risanje dvodimenzionalnih grafov, t. j. upodabljanje funkcij <math>z(x,y)</math>.

Naj bo funkcija podana s tabelo <tt>z[i,j]</tt>, <tt>i = 1, ..., nx</tt>, <tt>j = 1, ..., ny</tt>, pri čemer sta koordinati točke <math>(i,j)</math> na zaslonu <tt>ix = x0 + i*dx</tt>, <tt>jy = y0 + j*dy</tt>. Na znani način poiščemo zgornjo in spodnjo mejo, <tt>zmin</tt> in <tt>zmax</tt>. Recimo, da bi v sliki radi ločili (približno) 10 razredov vrednosti spremenljivke <math>z</math>. Najpreprostejša upodobitev je kar s senčenjem ali z barvo.

<tt>
 function barva(i, j: integer): integer;
 var zz: real;
 begin
   zz := 0.25*(z[i,j] + z[i+1,j] + z[i,j+1] + z[i+1,j+1]);
   barva := trunc(10.0*(zz - zmin)/(zmax - zmin))
 end;
</tt>

Narišemo pa:

<tt>
 for i := 1 to nx-1 do
   for j := 1 to ny-1 do
     rectangle(x0 + i*dx, y0 + j*dy, x0 + i*dx + dx, y0 + j*dy + dy, barva(i, j))
</tt>

Lepša je slika z izohipsami. Spet si bomo izbrali ~10 vrednosti za izohipse med <tt>zmin</tt> in <tt>zmax</tt>. Izohipso pri vrednosti <tt>z0</tt> narišemo tako, da pregledamo vse pravokotniške celice po vrsti, kot tečeta indeksa <tt>i</tt> in <tt>j</tt> njihovega spodnjega levega oglišča, ali vsebujejo vrednosti <tt>z0</tt>. Najprej pregledamo, kje na vsakem od robov celice je vrednost <tt>z0</tt>. Če vrednosti ni na tem robu, vrne funkcija vrednost <tt>-1</tt>, sicer pa najde točko na robu z linearno interpolacijo.

<tt>
 function nicla(i1, j1, i2, j2: integer; z0: real): real
 var z1, z2: real;
 begin
    z1 := z[i1,j1];
    z2 := z[i2,j2];
    if (z0 - z1)*(z0 - z2) < 0 then nicla := (z0 - z1)/(z2 - z1) else nicla := -1
 end;
</tt>

Zdaj pogledamo po vrsti robove celice in si beležimo, na kolikih od njih je vrednost z0. Možno je, da dobimo 2 točki ali 4. Ustrezno temu potegnemo en ali dva odseka izohipse.

<tt>
 procedure celica(i,j: integer; z0: real)
 var m: integer; p: real; xx, yy: array[1..4] of integer;
 begin
   m := 1; 
   p := nicla(i, j, i+1, j, z0);
   if p > 0 then begin 
     xx[m] := x0 + (i + p)*dx;
     yy[m] := y0 + j*dy; 
     inc(m)
   end;
   p := nicla(i+1, j, i+1, j+1, z0);
   if p > 0 then begin 
     xx[m] := x0 + (i + 1)*dx;
     yy[m] := y0 + (j + p)*dy;
     inc(m)
   end;
   p := nicla(i+1, j+1, i, j+1, z0);
   if p>0 then begin 
     xx[m] := x0 + (i + 1 - p)*dx;
     yy[m] := y0 + (j + 1)*dy;
     inc(m)
   end;
   p := nicla(i, j+1, i, j, z0);
   if p > 0 then begin  
     xx[m] := x0 + i*dx;
     yy[m] := y0 + (j + 1 - p)*dy;
   end;
   if m > 1 then plotline(xx[1], yy[1], xx[2], yy[2]);
   if m = 4 then plotline(xx[3], yy[3], xx[4], yy[4])
 end;
</tt>

V glavnem programu napravimo zanko po vseh vrednostih <tt>z0</tt>, za katere želimo narisati izohipse, in seveda dvojno zanko po vseh celicah:

<tt>
 for i := 1 to nx-1 do
   for j := 1 to ny-1 do
     celica(i, j, z0);
</tt>

Notranje izohipse bodo nekoliko oglate, zunanje pa dokaj dobre. Nekaj malega bi lahko izboljšali z uporabo odsekov višjega reda (n.pr. paraboličnih lokov), predvsem pa pomaga, če je tabela funkcijskih vrednosti <tt>z[i,j]</tt> dovolj gosta.

V '''Mathematici''' narišemo sliko z barvami ali izohipsami z ukazoma <tt>ListDensityPlot</tt> in <tt>ListContourPlot</tt>. Tabelo z 10 x 10 vrednostmi na intervalu <math>0<x<1</math> in <math>1<y<3</math> v '''gnuplotu''' prikažemo z desetimi barvami tako: 
<tt>
 set view map
 unset surface
 set pm3d
 set palette maxcolors 10
 splot "podatki.txt" matrix using ($1/9):(1+2*$2/9):3
</tt>
z desetimi izohipsami pa tako:
<tt>
 set view map
 unset surface
 set contour
 set cntrparam levels 10
 splot "podatki.txt" matrix using ($1/9):(1+2*$2/9):3 with lines
</tt>

== Naloge: ==

# Za 28 absorpcijskih spektrov robu K železa v datoteki "Fe_rob_0_27.xmu" nariši dvodimenzionalni graf, kjer je ena os energija fotona, druga pa čas v teku polnjenja in praznenja. Prikaži z barvo in izohipsami.
# Prikaži temperaturno polje v prečnem prerezu dimnika, kjer je temperatura vročih plinov 200 <sup>0</sup>C, na zunanji steni pa je 0 <sup>0</sup>C, iz podatkov v datoteki "Dimnik.dat". V datoteki je območje temperatur normirano na interval [0, 1], podane so v mreži 24 x 24 točk. Napravi grafa z barvno lestvico in z risanjem izoterm.
# Napravi graf izoterm <math>T(p,V)</math> za Van der Waalsov plin z enačbo stanja <math>(p+a/V^2)(V-b)=RT</math>, ki jo najprej predelamo v brezdimenzijsko obliko, tako da vse tri spremenljivke <math>p</math>, <math>V</math> in <math>T</math> normiramo na njihove vrednosti v kritični točki in se enačba v novih spremenljivkah <math>\Pi</math>, <math>\Phi</math> in <math>\Theta</math> glasi <math>(\Pi+3/\Phi^2)(3\Phi-1)=8\Theta</math>. Izberi primerno območje za spremenljivki <math>\Pi</math> in <math>\Phi</math> in si pripravi tabelo funkcije za risanje.