Cvičení č. 13: Vizualizace ve třech dimenzích

Poslední cvičení je věnováno přípravě a vizualizaci dat ve třech dimenzích. Půjde o banán, srdce, sférické harmoniky a atomové orbitaly. Pro kreslení užijeme Matplotlib, gnuplot, Octave a ParaView.

Parametrické funkce

Nejprve se budeme zabývat situací, kdy je objekt (funkce) ve třech dimenzích (3D) popsán parametricky. Máme na mysli zápis

\[x = x(u,v,w),\ y = y(u,v,w),\ z = z(u,v,w),\]

kde \(u,v,w\) jsou nezávislé parametry a \(x,y,z\) kartézské souřadnice bodů náležejících k 3D objektu. Při dvou nezávislých parametrech půjde o 2D plochu, při jediném o 1D křivku. V PDF souboru Plochy – potrava pro oko jsou shrnuty parametrické funkce popisující mj. banán, vajíčko, bagel, šneka či kapříka. Sem převezmeme banán:

\[x=a\sin\theta\cos\phi,\ y=b\sin\theta\sin\phi,\ z=b(\cos\theta+c\sin^2\theta\cos^2\phi),\]

kde \(0\le\theta\le\pi\) a \(0\le\phi\le2\pi\) jsou nezávislé parametry a ostatní fixujeme na hodnotách \(a=1\), \(b=1/5\) a \(c=2\).

Pomocí vnořených cyklů systematicky projdeme požadované intervaly nezávislých parametrů, vytvoříme tak síť uzlových bodů a zapíšeme trojice jejich kartézských souřadnic. K uzlům bychom mohli přidávat další (skalární nebo vektorová) data a ty pak – kromě samotného tvaru objektu – také vizualizovat (např. barevně nebo šipkami). Pokud by se tato data dala vyjádřit jen ze souřadnic, mohli bychom jejich výpočet ponechat až na vizualizační program. My budeme chtít obarvit banán realisticky na žluto a na obou koncích do černa. Uhodneme, že pro konce platí \(\sqrt{x^2+y^2+z^2}\gtrsim a\), a ke každému uzlu přidáme velikost jeho průvodiče, kterou navážeme na vhodnou barevnou paletu, takovou, aby pro hodnoty menší než \(a\) byla žlutá a pro hodnoty větší než \(a\) přešla do černé. Pro gnuplot budeme generovat datový soubor jen se souřadnicemi, velikost průvodiče dopočteme v gnuplotu. Pro ParaView zapíšeme VTK soubor (ve variantě STRUCTURED_GRID pro křivočarou síť) s hlavičkou, sekcí souřadnic a sekcí dat. Komentované ukázky, jak mají VTK soubory vypadat, obsahuje PDF soubor Poprvé s ParaView, a vykoukat to snad lze i z následujících zdrojových textů.

# Python: příprava dat pro vizualizaci banánu v gnuplotu, ParaView a Matplotlibu
import numpy as np

# Vytvoření polí s kartézskými souřadnicemi a daty
def createBanana(nth,nph):
  a=1; b=0.2; c=2                 # parametry banánu
  nr=0                            # 0 pro jednu plochu
  thmin=0; thmax=np.pi            # meze souřadnice theta (poledník): theta_min a theta_max
  phmin=0; phmax=np.pi*2          # meze souřadnice phi (rovnoběžka): phi_min a phi_max
  x=np.empty((nr+1,nth+1,nph+1))  # pole kartézských souřadnic
  y=x.copy()
  z=x.copy()
  s=x.copy()                      # skalární data v uzlech sítě
  ir=nr                           # jedna plocha
  for ith in range(nth+1):        # v NumPy měníme indexy nejrychleji zprava
    th=thmin+ith*(thmax-thmin)/nth
    for iph in range(nph+1):
      ph=phmin+iph*(phmax-phmin)/nph
      x[ir,ith,iph]=a*np.sin(th)*np.cos(ph)
      y[ir,ith,iph]=b*np.sin(th)*np.sin(ph)
      z[ir,ith,iph]=b*(np.cos(th)+c*(np.sin(th)*np.cos(ph))**2)
  s=np.sqrt(x*x+y*y+z*z)
  return (x,y,z,s)

# Vytvoření PNG-souboru v Matplotlibu (plocha v 3D)
def saveFig(x,y,z,s,file):
  import matplotlib.pyplot as plt
  import matplotlib.cm as cm
  import matplotlib.colors as colors
  with open(file,'w') as f:
    fig=plt.figure(layout='tight')
    ax=plt.axes(projection='3d')
    cmap=colors.ListedColormap(['yellow','yellow','yellow','black'])
    # ax.plot_surface(x[0,:,:],y[0,:,:],z[0,:,:],cmap=cm.inferno)  # z-coordinate by predefined colormap
    # ax.plot_surface(x[0,:,:],y[0,:,:],z[0,:,:],cmap=cmap)        # z-coordinate by defined cmap
    ax.plot_surface(x[0,:,:],y[0,:,:],z[0,:,:],facecolors=cmap(s[0,:,:])) # s-array by defined cmap
    ax.set_aspect('equal')
    ax.view_init(elev=10,azim=120)
    plt.savefig(file)
    plt.show()

# Vytvoření DAT-souboru pro gnuplot (plocha v 3D)
def saveDat(x,y,z,file):
  (nr,nth,nph)=np.shape(x)   # počty uzlů
  fmt='8.4f'                 # formát výstupu
  with open(file,'w') as f:
    for ir in range(nr):
      for ith in range(nth):
        for iph in range(nph):
          print(f'{x[ir,ith,iph]:{fmt}} {y[ir,ith,iph]:{fmt}} {z[ir,ith,iph]:{fmt}}',file=f)

# Vytvoření VTK-souboru pro ParaView (plocha v 3D)
# varianta STRUCTURED_GRID pro křivočarou síť/curvilinear structured grid
def saveVtkSG(x,y,z,s,file):
  (nr,nth,nph)=np.shape(x)   # počty uzlů
  fmt='8.4f'                 # formát výstupu
  with open(file,'w') as f:
    print('# vtk DataFile Version 3.0',file=f)
    print('vtk output',file=f)
    print('ASCII',file=f)
    print('DATASET STRUCTURED_GRID',file=f)
    print('DIMENSIONS',nr,nth,nph,file=f)
    print('POINTS',nr*nth*nph,'float',file=f)
    for iph in range(nph):   # pro ParaView měníme indexy nejrychleji zleva
      for ith in range(nth):
        for ir in range(nr):
          print(f'{x[ir,ith,iph]:{fmt}} {y[ir,ith,iph]:{fmt}} {z[ir,ith,iph]:{fmt}}',file=f)
    print('POINT_DATA',nr*nth*nph,file=f)
    print('SCALARS MyData float 1',file=f)
    print('LOOKUP_TABLE default',file=f)
    for iph in range(nph):
      for ith in range(nth):
        for ir in range(nr):
          print(f'{s[ir,ith,iph]:{fmt}}',file=f)

nth=90   # počet kroků sférické sítě ve směru theta (poledník)
nph=180  # počet kroků sférické sítě ve směru phi (rovnoběžka)
(x,y,z,s)=createBanana(nth,nph)
if True: soubor='banan.dat'; saveDat(x,y,z,soubor);     print('Soubor',soubor,'uložen.')
if True: soubor='banan.vtk'; saveVtkSG(x,y,z,s,soubor); print('Soubor',soubor,'uložen.')
if True: soubor='banan.png'; saveFig(x,y,z,s,soubor);   print('Soubor',soubor,'uložen.')

Datový soubor banan.dat pro gnuplot obsahuje pouze tři sloupce se souřadnicemi síťových bodů. Velikost průvodiče v bodech sítě vytvoříme až v gnuplotu příkazem splot jako virtuální čtvrtý sloupec, klauzulí palette zohledněný pro obarvení jednotlivých uzlů. Žlutočernou barevnou paletu si připravíme příkazem set palette defined.

# gnuplot: banán
unset border
unset xtics
unset ytics
unset ztics
unset colorbox
unset key
set palette defined (0 'yellow', 0.8 'yellow', 1.1 'black')
set view 75,135
splot [-1:1][-1:1][-1:1] 'banan.dat' using 1:2:3:(sqrt(column(1)**2+column(2)**2+column(3)**2)) palette with lines
set terminal pngc size 2500,1500 crop
set output 'banan-gnuplot.png'
replot
set output
unset terminal

Banán v gnuplotu…

_images/banan-gnuplot.png

S přichystaným souborem banan.vtk lze v ParaView už jen klikat myší: File–Open, načíst soubor banan.vtk, Apply, pootočit banánem při stisknutém levém tlačítku myši (a zkusit, co dělá pohyb myši při stisku ostatních tlačítek), nahoře ikonkami Toggle Color Legend Visibility a Show Orientation Axes vypnout barevnou škálu a osový kříž, ikonkou Load a color palette nastavit White Background, vlevo v části Coloring stisknout Edit pro Color Map Editor, v obrázku Mapping Data dvojkliknout na levém kolečku ve vodorovném pruhu barevné škály a změnit barvu na žlutou, střední kolečko posunout někam pod 1,0 a změnit také na žlutou, pravému kolečku nastavit černou, pak ještě vlevo v Lighting pro parádu přidat na Specular a konečně ve File–Save Screenshot nastavit velikost výstupního PNG obrázku a uložit ho.

Banán v ParaView…

_images/banan-paraview.png

Implicitní funkce

Dále se zaměříme na případ, kdy je plocha popsána implicitně rovnicí

\[f(x,y,z)=0.\]

Ze souboru Plochy – potrava pro oko (a odsud) nyní vybereme rovnici popisující 3D srdce:

\[\left(ax^2+by^2+cz^2-1\right)^3-dx^2z^3-ey^2z^3=0,\]

kde sada konstant \(a=1\), \(b=9/4\), \(c=1\), \(d=1\), \(e=9/80\) definuje spíše široké srdce, zatímco sada \(a=2\), \(b=2\), \(c=1\), \(d=1\), \(e=1/10\) popisuje úzký exemplář.

Běžnou cestou k vizualizaci implicitně zadané plochy je tabelování funkce \(f(x,y,z)\) na jednoduché 3D oblasti (např. na kvádru) a zobrazení její (nulové) izoplochy. Tuto schopnost nemá gnuplot, má ji však Octave (funkce isosurface) i ParaView (modul Contour).

Skript pro Octave (verze 6, s verzí 7 je něco špatně), který na kartézské síti [x,y,z] funkci \(f(x,y,z)\) tabeluje do pole s a rovnou i kreslí její izoplochu, isosurface(x,y,z,s,0), následuje. Je použitelný i v MATLABu, který – na rozdíl od Octavu – umí zacházet i s odlesky, což se snadno vyzkouší nahrazením aktivního řádku patch zakomentovanou alternativou.

% Octave: 3D srdce
a=1; b=9/4; c=1; d=1; e=9/80;
xmin=-1.2; xmax=1.2; ymin=-.75; ymax=.75; zmin=-1.2; zmax=1.5;
nx=100; ny=100; nz=300;
xvec=linspace(xmin,xmax,nx+1); yvec=linspace(ymin,ymax,ny+1); zvec=linspace(zmin,zmax,nz+1);
[x,y,z]=ndgrid(xvec,yvec,zvec);
s=(a*x.^2+b*y.^2+c*z.^2-1).^3-d*x.^2.*z.^3-e*y.^2.*z.^3;
set(gcf,'Color','white'); clf; axis off; view([-1 5 0.5]); % bílé pozadí, bez os, pozorovací bod
[f,v]=isosurface(x,y,z,s,0);
patch('Faces',f,'Vertices',v,'FaceColor','red','EdgeColor',[.5 .0 .0]);
% MATLAB (with lighting)
% patch('Faces',f,'Vertices',v,'FaceColor','red','EdgeColor','none'); camlight('right'); lighting gouraud;

Srdce v Octavu (3D dojem vytvářejí hrany elementů tvořících izoplochu) a MATLABu (3D dojem díky stínování)…

_images/srdce-octave.png _images/srdce-matlab.png

Datový soubor pro ParaView s funkcí \(f(x,y,z)\) tabelovanou ve VTK formátu (varianta STRUCTURED_POINTS pro rovnoměrnou síť, viz Poprvé s ParaView) připravíme pomocí následujícího programu. Tabelujeme pole s na síti x,y,z, stejně jako v předešlém skriptu. Přípravu izoplochy necháváme na ParaView.

# Python: příprava dat pro vizualizaci srdce v ParaView
import numpy as np

# Vytvoření polí s kartézskými souřadnicemi a daty
def createHearth(nx,ny,nz):
  a=1; b=9/4; c=1; d=1; e=9/80    # parametry širokého srdce
  # a=2; b=2; c=1; d=1; e=1/10    # parametry úzkého srdce
  xmin=-1.2; xmax=1.2             # meze ve směru x
  ymin=-0.75;ymax=0.75            # meze ve směru y
  zmin=-1.2; zmax=1.5             # meze ve směru z
  s=np.empty((nx+1,ny+1,nz+1))    # skalární data v uzlech sítě
  for ix in range(nx+1):          # v NumPy měníme indexy nejrychleji zprava
    x=xmin+ix*(xmax-xmin)/nx
    for iy in range(ny+1):
      y=ymin+iy*(ymax-ymin)/ny
      for iz in range(nz+1):
        z=zmin+iz*(zmax-zmin)/nz
        s[ix,iy,iz]=(a*x**2+b*y**2+c*z**2-1)**3-d*x**2*z**3-e*y**2*z**3
  return ((xmin,xmax,ymin,ymax,zmin,zmax),s)

# Vytvoření VTK-souboru pro ParaView (3D data v kvádru)
# varianta STRUCTURED_POINTS pro rovnoměrnou kartézskou síť/uniform rectilinear grid
def saveVtkSP(bounds,s,file):
  (xmin,xmax,ymin,ymax,zmin,zmax)=bounds
  (nx,ny,nz)=np.shape(s)       # počty uzlů
  fmt='15.10e'                 # formát výstupu
  with open(file,'w') as f:
    print('# vtk DataFile Version 3.0',file=f)
    print('vtk output',file=f)
    print('ASCII',file=f)
    print('DATASET STRUCTURED_POINTS',file=f)
    print('DIMENSIONS',nx,ny,nz,file=f)
    print('ORIGIN',xmin,ymin,zmin,file=f)
    print('SPACING',(xmax-xmin)/(nx-1),(ymax-ymin)/(ny-1),(zmax-zmin)/(nz-1),file=f)
    print('POINT_DATA',nx*ny*nz,file=f)
    print('SCALARS MyData float 1',file=f)
    print('LOOKUP_TABLE default',file=f)
    for iz in range(nz):     # pro ParaView měníme indexy nejrychleji zleva
      for iy in range(ny):
        for ix in range(nx):
          print(f'{s[ix,iy,iz]:{fmt}}',file=f)

nx=100  # počet uzlů kartézské sítě ve směru x
ny=100  # počet uzlů kartézské sítě ve směru y
nz=100  # počet uzlů kartézské sítě ve směru z
(bounds,s)=createHearth(nx,ny,nz)
if True: soubor='srdce.vtk'; saveVtkSP(bounds,s,soubor); print('Soubor',soubor,'uložen.')

Izoplochy v ParaView: File–Open, načíst soubor srdce.vtk, Apply, vlevo nahoře kliknout na ikonu Contour (nebo ekvivalentně volit Filters–Common–Contour), vlevo v Isosurfaces změnit hodnotu na 0, Apply, nahoře ikonou Set view direction to +Y nebo myší nastavit vhodnější úhel pohledu, vlevo v Coloring–Edit vybrat červenou barvu, odstranit bounding box kliknutím na ikoně oka vlevo u souboru srdce.vtk a stejně jako u banánu vybrat barvu pozadí, nastavit odrazivost a uložit obrázek.

Srdce v ParaView…

_images/srdce-paraview-1.png _images/srdce-paraview-2.png

Můžeme ovšem ještě pokračovat: View–Animation View, v záhlaví okna nastavit hodnoty Start Time a End Time třeba na -10 a 10, No. Frames pak řekněme 101, dále modré plus u Contour1–Isosurfaces, tam pak dvojklik pro Animation Keyframes, nastavit shodné meze pro čas (-10, 10) a hodnoty Value od -.05 do +.05, OK a zelená šipka pro Play. Po File–Save Animation získáme toto čili animaci izoploch blízkých nulové izoploše.

Sférické harmonické funkce

Se sférickými harmonikami \(Y_l^m(\theta,\phi)\) se potkáváme v situacích, kde hraje roli sférická symetrie (mj. v jaderné fyzice, kvantové fyzice, astronomii nebo geofyzice). Sférické harmoniky jsou ortonormální na povrchu koule, funkce definované na kouli lze reprezentovat jako řadu sférických harmonik, funkce \(r^lY_l^m(\theta,\phi)\) a \(r^{-l-1}Y_l^m(\theta,\phi)\) vyhovují Laplaceově rovnici v 3D, \(\Delta f(r,\theta,\phi)=0\). Ať tak či onak, jejich matematická definice zní:

\[Y_l^m(\theta,\phi)=N_l^m P_l^m(\cos\theta)e^{im\phi},\]

kde \(\theta,\ \phi\) jsou úhlové sférické souřadnice, normovací konstanty \(N_l^m=(-1)^m\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}\) a \(P_l^m(\cos\theta)\) značí přidružené Legendrovy funkce. Ty jsou zavedeny (a ortogonální) na intervalu \(\langle-1,1\rangle\) a numericky se vyčíslují pomocí rekurentního vztahu:

\[P_{k+2}^m(x)=\frac{(2k+3)xP_{k+1}^m(x)-(k+m+1)P_k^m(x)}{k-m+2},\ k=m,m+1,\dots,\]

který se startuje pomocí vzorců \(P_m^m(x)=(-1)^m\frac{(2m)!}{2^m m!}(1-x^2)^{\frac m2}\) a \(P_{m+1}^m(x)=(2m+1)xP_m^m(x)\). Tyto tři vztahy zobecňují rekurentní vzorec pro Legendrovy polynomy \(P_l(x)\equiv P_l^0(x)\), s nimiž jsme se setkali na přednášce o numerickém integrování: \(P_{k+2}(x)=\frac{(2k+3)xP_{k+1}(x)-(k+1)P_k(x)}{k+2},\ k=0,1,\dots,\) se startovacími hodnotami \(P_0(x)=1,\ P_1(x)=x\).

# Python: příprava dat pro vizualizaci sférických harmonických funkcí
import numpy as np
import scipy

# Vytvoření polí s kartézskými souřadnicemi a daty
def createSphHarm(nth,nph):
  n=3; m=3                        # stupeň a řád sférické harmonické funkce
  nr=0                            # 0 pro jednu plochu
  thmin=0; thmax=np.pi            # meze souřadnice theta (poledník): theta_min a theta_max
  phmin=0; phmax=np.pi*2          # meze souřadnice phi (rovnoběžka): phi_min a phi_max
  x=np.empty((nr+1,nth+1,nph+1))  # pole kartézských souřadnic
  y=x.copy()
  z=x.copy()
  s=x.copy()                      # skalární data v uzlech sítě
  ir=nr                           # jedna plocha
  for ith in range(nth+1):        # v NumPy měníme indexy nejrychleji zprava
    th=thmin+ith*(thmax-thmin)/nth
    for iph in range(nph+1):
      ph=phmin+iph*(phmax-phmin)/nph
      s[ir,ith,iph]=scipy.special.sph_harm(m,n,ph,th).real
#     ss=1
      ss=abs(s[ir,ith,iph])       # faktor pro lokální deformaci kulové plochy
#     ss=s[ir,ith,iph]**2
      x[ir,ith,iph]=ss*np.sin(th)*np.cos(ph)
      y[ir,ith,iph]=ss*np.sin(th)*np.sin(ph)
      z[ir,ith,iph]=ss*np.cos(th)
  return (x,y,z,s)

# Vytvoření PNG-souboru v Matplotlibu (plocha v 3D)
def saveFig(x,y,z,s,file):
  import matplotlib.pyplot as plt
  import matplotlib.cm as cm
  import matplotlib.colors as colors
  with open(file,'w') as f:
    fig=plt.figure(layout='tight')
    ax=plt.axes(projection='3d')
    cmap=colors.ListedColormap(['yellow','yellow','yellow','black'])
    ax.plot_surface(x[0,:,:],y[0,:,:],z[0,:,:],cmap=cm.inferno)  # z-coordinate by predefined colormap
    # ax.plot_surface(x[0,:,:],y[0,:,:],z[0,:,:],cmap=cmap)        # z-coordinate by defined cmap
    # ax.plot_surface(x[0,:,:],y[0,:,:],z[0,:,:],facecolors=cmap(s[0,:,:])) # s-array by defined cmap
    ax.set_aspect('equal')
    ax.view_init(elev=10,azim=120)
    plt.savefig(file)
    plt.show()

# Vytvoření VTK-souboru pro ParaView (plocha v 3D)
# varianta STRUCTURED_GRID pro křivočarou síť/curvilinear structured grid
def saveVtkSG(x,y,z,s,file):
  (nr,nth,nph)=np.shape(x)   # počty uzlů
  fmt='8.4f'                 # formát výstupu
  with open(file,'w') as f:
    print('# vtk DataFile Version 3.0',file=f)
    print('vtk output',file=f)
    print('ASCII',file=f)
    print('DATASET STRUCTURED_GRID',file=f)
    print('DIMENSIONS',nr,nth,nph,file=f)
    print('POINTS',nr*nth*nph,'float',file=f)
    for iph in range(nph):   # pro ParaView měníme indexy nejrychleji zleva
      for ith in range(nth):
        for ir in range(nr):
          print(f'{x[ir,ith,iph]:{fmt}} {y[ir,ith,iph]:{fmt}} {z[ir,ith,iph]:{fmt}}',file=f)
    print('POINT_DATA',nr*nth*nph,file=f)
    print('SCALARS MyData float 1',file=f)
    print('LOOKUP_TABLE default',file=f)
    for iph in range(nph):
      for ith in range(nth):
        for ir in range(nr):
          print(f'{s[ir,ith,iph]:{fmt}}',file=f)

nth=90   # počet kroků sférické sítě ve směru theta (poledník)
nph=180  # počet kroků sférické sítě ve směru phi (rovnoběžka)
(x,y,z,s)=createSphHarm(nth,nph)
if True: soubor='sh.vtk'; saveVtkSG(x,y,z,s,soubor); print('Soubor',soubor,'uložen.')
if True: soubor='sh.png'; saveFig(x,y,z,s,soubor);   print('Soubor',soubor,'uložen.')

Sférické harmoniky můžeme vizualizovat poctivě na povrchu koule vynesením barvy přidružené hodnotě jejich reálné části, tedy barevnou reprezentací skalárního pole \(\operatorname{Re}(Y_l^m(\theta,\phi))\) na ploše popsané parametricky:

\[x = \sin\theta\cos\phi,\ y = \sin\theta\sin\phi,\ z = \cos\theta.\]

Zonální sférické harmoniky \(Y_l^0(\theta,\phi),\ l=0,1,2,3,\) na kulové ploše…

\[\textstyle Y_0^0=N_0^0,\ Y_1^0=N_1^0\cos\theta,\ Y_2^0=N_2^0\left(\frac32\cos^2\theta-\frac12\right),\ Y_3^0=N_3^0\left(\frac52\cos^3\theta-\frac32\cos\theta\right)\]
_images/y00-bw.png _images/y10-bw.png _images/y20-bw.png _images/y30-bw.png

Mnozí to však považují za nudné a obarvenou kulovou plochu, \(r=1,\) ještě přetvarují podle předpisu \(r=|\operatorname{Re}(Y_l^m(\theta,\phi))|\) nebo \(r=|\operatorname{Re}(Y_l^m(\theta,\phi))|^2\):

\[x = |\operatorname{Re}(Y_l^m(\theta,\phi))|\sin\theta\cos\phi,\ y = |\operatorname{Re}(Y_l^m(\theta,\phi))|\sin\theta\sin\phi,\ z = |\operatorname{Re}(Y_l^m(\theta,\phi))|\cos\theta.\]

Zonální sférické harmoniky \(Y_l^0(\theta,\phi),\ l=0,1,2,3,\) v ParaView (osa z míří vzhůru mírně napravo)…

\[\textstyle Y_0^0=N_0^0,\ Y_1^0=N_1^0\cos\theta,\ Y_2^0=N_2^0\left(\frac32\cos^2\theta-\frac12\right),\ Y_3^0=N_3^0\left(\frac52\cos^3\theta-\frac32\cos\theta\right)\]
_images/y00.png _images/y10.png _images/y20.png _images/y30.png

Reálné sektorové sférické harmoniky \(\operatorname{Re}(Y_l^l(\theta,\phi)),\ l=0,1,2,3\) (osa z míří doprava mírně vpřed)…

\[\textstyle Y_0^0=N_0^0,\ Y_1^1=N_1^1(-\sin\theta)e^{i\phi},\ Y_2^2=N_2^2(3\sin^2\theta)e^{2i\phi},\ Y_3^3=N_3^3(-15\sin^3\theta)e^{3i\phi}\]
_images/y00.png _images/y11.png _images/y22.png _images/y33.png

…a jejich kombinace – proč ji nenazvat kapřík:

_images/yll-kaprik1.png _images/yll-kaprik2.png _images/yll-kaprik3.png

Příprava dat i práce s ParaView je tu shodná jako při malování banánu neboli parametricky zadané funkce (o dva odstavce výše). Vlastnoručně zhotovená data malovaná v ParaView bychom mohli směle porovnávat s obrázkem na sférických harmonik na Wikipedii.

Atomové orbitaly

Atomový orbital znázorňuje pravděpodobnost výskytu elektronu v okolí atomového jádra. Ta se odvozuje od 3D vlnové funkce \(\psi_{nlm}(r,\theta,\phi)\), která má pro vodíkový atom následující tvar (odvození např. zde):

\[\psi_{nlm}(r,\theta,\phi)=\sqrt{\left(\frac{2}{na}\right)^3 \frac{(n-l-1)!}{2n(n+l)!}} \ e^{-\frac{r}{na}} \ \left(\frac{2r}{na}\right)^l \ L_{n-l-1}^{2l+1}\left(\frac{2r}{na}\right) \ Y_l^m(\theta,\phi).\]

Indexy \(n,l,m\) se nazývají hlavní, vedlejší a magnetické kvantové číslo, konstanta \(a\) = 5,292.10-11 m je Bohrův poloměr, \(L_k^\alpha(x)\) značí zobecněné Laguerrovy polynomy a \(Y_l^m(\theta,\phi)\) sférické harmonické funkce, diskutované o odstavec výše. Zobecněné Laguerrovy polynomy, ortogonální s vahou \(x^\alpha e^{-x}\) na intervalu \(\langle0,\infty)\), vyčíslíme pomocí rekurentního vztahu:

\[L_{k+2}^\alpha(x)=\frac{(2k+3+\alpha-x)L_{k+1}^\alpha(x)-(k+1+\alpha)L_k^\alpha(x)}{k+2},\ k=0,1,\dots,\]

startovaného z hodnot \(L_0^\alpha(x)=1\) a \(L_1^\alpha(x)=1+\alpha-x\).

Pro zobrazení atomového orbitalu pro kvantová čísla \(n,l,m\) je třeba vyčíslit na vhodně zvolené oblasti vlnovou funkci \(\psi_{nlm}(r,\theta,\phi)\), přejít k její reálné reprezentaci (vhodnou lineární kombinací komplexních vlnových funkcí nebo jednoduše omezením se na \(\operatorname{Re}(\psi_{nlm})\)), spočítat pravděpodobnost \(P(r,\theta,\phi)\) jako kvadrát reálné vlnové funkce a vybrat její izoplochu, \(P(r,\theta,\phi)=\mathrm{const},\) která ohraničuje oblast s dostatečně vysokou pravděpodobností výskytu elektronu nebo prostě jen dobře vypadá. (Vytvořit izoplochu v ParaView jsme si už zkusili na příkladu 3D srdce.) Na Wikipedii vizualizují orbitaly pomocí samotné vlnové funkce (nikoliv pravděpodobnosti), tak uděláme totéž, abychom mohli srovnávat.

Orbitaly vodíkového atomu: izoplochy \(\psi_{210},\ \psi_{310},\ \psi_{320},\ \psi_{321}\) v ParaView…

_images/orb-210.png _images/orb-310.png _images/orb-320.png _images/orb-321.png

Orbitaly \(\psi_{410},\ \psi_{420},\ \psi_{421},\ \psi_{432}\)

_images/orb-410.png _images/orb-420.png _images/orb-421.png _images/orb-432.png

Orbitaly \(\psi_{510},\ \psi_{520},\ \psi_{530},\ \psi_{540}\)

_images/orb-510.png _images/orb-520.png _images/orb-530.png _images/orb-540.png

Program pro přípravu atomových orbitalů by se podobal programu pro přípravu srdce neboli izoplochy implicitní funkce (o dva odstavce výše). Jistě by zhodnotil funkce Leg a Lag z výše uvedených ukázek. A jsme u posledního domácího úkolu: sestavit program pro tabelování vlnové funkce atomu vodíku na vhodně voleném kvádru a pomocí ParaView získat několik vizuálně kvalitních orbitalů.