Využití IDL pro zobrazení výsledků numerických výpočtů, program shell.pro

Richard Wunsch, 1. ročník, PgS fyzika

Program se ve skutečnosti skládá z několika procedur, které slouží k vytvoření obrázků v publikační kvalitě k článku Wunsch R., Palous J.: "Gravitational instability of expanding shells: Solution with nonlinear terms", Astronomy & Astrophysics, 2001, submitted.

HI obálky jsou struktury v mezihvězdné hmotě v galaxiích. Vznikají při uvolnění velkého množství energie - řádově 10⁵¹-10⁵⁴ ergů. Mechanismem může být např. činnost OB asociace (exploze supernov a hvězdný vítr), dopad trpasličí galaxie nebo vysokorychlostního plynného oblaku na galaktický disk, nebo událost, která vede ke vzniku gamma záblesku. Vznikne tak jakási bublina o poloměru 10-1000 pc vyplněná řídkým horkým plynem, téměř veškerá hmota je soustředěna na jejím povrchu. Obálka se rozpíná a akreuje na sebe další hmotu a může se stát gravitačně nestabilní a začít fragmentovat.

Článek se zabývá nelineární analýzou gravitačních nestabilit těchto obálek. Výsledkem úprav hydrodynamických rovnic a Poissonovy rovnice je soustava tří diferenciálních rovnic pro tři komplexní veličiny (amplitudy), které popisují časový vývoj tří nejméně stabilních modů povrchové hustoty, které mají všechny stejnou velikost a jsou vzájemně skloněny o 60 stupňů. Pro povrchovou hustotu tedy máme:

sigma(Theta)   = sig(xi) * exp(i*(eta,Theta)) 
               + sig(xi₊) * exp(i*(eta₊,Theta))    
               + sig(xi_-) * exp(i*(eta_-,Theta))

kde sig() je funkce, která přepočítává bezrozměrné komplexní amplitudy xi, xi₊ a xi_- na amplitudy hustoty, eta, eta₊ a eta_- jsou vlnové vektory (vzájemně skloněné o 60 stupňů) odpovídajících modů a Theta = (Theta_x, Theta_y) jsou souřadnice na povrchu obálky.

Výsledkem numerického řešení diferenciálních rovnic je soubor, který obsahuje časový vývoj amplitud jednotlivých modů (první sloupec čas, dalších šest sloupců hodnoty amplitud).

POVRCHOVÁ HUSTOTA

První část programu - procedura sig_surf spočítá a zobrazí časový vývoj rozložení povrchové hustoty. Skládá se zhruba ze tří částí:

Volání procedury read_data, která načte soubory s důležitými daty - především soubor xi - hodnoty bezrozměrných amplitud tří interagujících modů, dále soubory physical a math, které obsahují další důležité údaje o expandující obálce (poloměr, expanzní rychlost apod.; používá je např. funkce sig()).
Volání procedury make_sig, která z amplitud xi vypočítá povrchovou hustotu. Výsledek uloží do třírozměrného pole s_field (čas + 2 prostorové souřadnice). Při výpočtu se volají funkce get_sigma_1, která vrátí povrchovou hustotu v daném čase a bodě, a funkce get_sigma_eta, get_sigma_epl, get_sigma_emi, které přepočítávají amplitudy xi na povrchovou hustotu - implementace funkce sig().
Vlastní vykreslení povrchové hustoty. K vykreslení je použita knihovní procedura surface v kombinaci s procedurou contour, která kreslí obrázek transformovaný tak, aby odpovídal rovině XY obrázku nakresleného procedurou surface - to zajišťuje nastavení klíčového parametru /t3d. Obrázky se ukládají do postscriptového souboru - každý časový krok jedna stránka. Nakonec se ještě vytvoří soubor fig1.eps s povrchovou hustotou v jednom konkrétním čase (50 Myr), který bude vložen do textu článku.

fig1.eps

RYCHLOSTNÍ POLE

Z amplitud xi, xi₊ a xi_- je možné vypočítat rychlostní pole toků hmoty na povrchu obálky:

velocity(Theta)   = vel(xi) * exp(i*(eta,Theta)) 
                  + vel(xi₊) * exp(i*(eta₊,Theta)) 
	          + vel(xi_-) * exp(i*(eta_-,Theta))

Veličiny mají podobný význam jako v případě povrchové hustoty, jen funkce vel() vrací 2D vektor.

Výpočet zajišťuje procedura vel_field, která podobně jako sig_surf nejdříve zavolá read_data, pak make_vel, která s pomocí funkcí get_v, get_v_eta, get_v_epl, get_v_emi vypočítá rychlostní pole v_field (tentokrát je čtyřrozměrné - čas, 2 souřadnice polohy, x nebo y složka rychlosti) a nakonec toto rychlostní pole vykreslí.

K vykreslení je tentokrát použita knihovní procedura velovect spolu s procedurou contour, která do obrázku přidá pro porovnání "vrstevnice" hustoty. Zde se objevil drobný problém se správným sesazením obrázků do sebe, protože procedura velovect, která se stará o nakreslení rychlostního pole, není zřejmě zcela standardní součástí IDL (alespoň verze 4) a nerozumí běžným grafickým klíčovým slovům jako např. [xy]range, [xy]style apod.

fig2.eps

FRAGMENTAČNÍ INTEGRÁL

Lineární teorie zavádí pro popis nestabilit tzv. rychlost růstu perturbací omega, která je obecně funkcí času a velikosti vlnového vektoru eta a má následující význam: mody povrchové hustoty s daným vlnovým číslem rostou v daném čase jako exp(omega*t):

sigma(eta, t) ~ exp(omega(eta,t)*t)

Pro expandující obálky se dá omega určit jako funkce poloměru R(t) a expanzní rychlosti V(t) a vlnového čísla eta.

Definuje se také tzv. fragmentační integrál:

kde t_b je začátek nestability. Fragmentační integrál tedy znamená něco ve smyslu "jak dlouho a jak hodně" byla obálka nestabilní. Jestliže tato bezrozměrná veličina dosáhne pro určitou vlnovou délku jisté rozumné hodnoty (např. 1), můžeme říci, že fragmenty této vlnové délky měly dostatek času se vyvinout. Jedním z cílů článku bylo prověřit toto poněkud volněji definované kritérium přesnější nelineární analýzou. Proto bylo potřeba zobrazit fragmentační integrál jako funkci eta a času.

Jeho výpočet a vykreslení provádí procedura i_frag. Nejdříve volá read_data, pak make_oe, která vypočítá z fyzikálních parametrů (R, V, apod) omega(eta, t) -> výsledek uloží do 2D pole omega_eta. Dále se volá procedura make_ifrag, která provede vlastní integraci s pomocí knihovní procedury IDL int_tabulated. Nakonec se pomocí procedury contour vykreslí výsledek.

fig3.eps

Závěrem něco o technických podmínkách

Celou práci jsem prováděl na pracovní stanici Sun s IDL verze 4 vzdáleně z Linuxového počítače. Proto byl kladen důraz na správné vykreslení obrázků do postscriptových souborů a cokoli se kreslí na obrazovku je vpodstatě jen z ladicích důvodů. Jednotlivé části programu vznikaly tak, jak přicházely a měnily se požadavky na obrázky do článku, takže nemá žádnou ucelenou koncepci a obsahuje spoustu provizorních řešení. Proto ho prosím berte trochu shovívavě, kdyby se přesto někdo chtěl na něco zeptat, může mě kontaktovat na adrese richard@saiph.ig.cas.cz.

R.W. , poslední úpravy 13.3.2001