8. Razvejitve

Iz Računalniška orodja v fiziki 2008 - 2009

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 (for i := 1 to n do).

Gibkost programiranja pa zagotavljajo predvsem razvejitveni ukazi, ko je nadaljnji potek računanja odvisen od izpolnitve kontrolnega pogoja (if pogoj then do ... else do ...). 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 z(x,y).

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

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;

Narišemo pa:

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))

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

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;

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.

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;

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

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

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 z[i,j] dovolj gosta.

V Mathematici narišemo sliko z barvami ali izohipsami z ukazoma ListDensityPlot in ListContourPlot. Tabelo z 10 x 10 vrednostmi na intervalu 0 < x < 1 in 1 < y < 3 v gnuplotu prikažemo z desetimi barvami tako:

set view map
unset surface
set pm3d
set palette maxcolors 10
splot "podatki.txt" matrix using ($1/9):(1+2*$2/9):3

z desetimi izohipsami pa tako:

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

Naloge:

1. 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.
2. Prikaži temperaturno polje v prečnem prerezu dimnika, kjer je temperatura vročih plinov 200 ⁰C, na zunanji steni pa je 0 ⁰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.
3. Napravi graf izoterm T(p,V) za Van der Waalsov plin z enačbo stanja (p + a / V²)(V − b) = RT, ki jo najprej predelamo v brezdimenzijsko obliko, tako da vse tri spremenljivke p, V in T normiramo na njihove vrednosti v kritični točki in se enačba v novih spremenljivkah Π, Φ in Θ glasi (Π + 3 / Φ²)(3Φ − 1) = 8Θ. Izberi primerno območje za spremenljivki Π in Φ in si pripravi tabelo funkcije za risanje.