Programovacie úlohy k predmetu „Metody zpracováni astrofyzikálních dat"
Mgr. Miroslav Jagelka 2020
Kapitola 1
Jednoduché programovacie úlohy
1. Určte tromi rôznymi spôsobmi inverznú maticu B k matici A.
2. Znázornite v grafe funkciu cos x v intervale [—27r,27r].
3. Vytvorte maticu M, v ktorej každý prvok bude obsahovať hodnotu dĺžky přepony zostrojenej nad pravouhlým trojuholníkom o odvěsnách s dĺžkami rovnými riadkovému a stĺpcovému indexu daného prvku matice. Matica nech je štvorcová o rozmere N.
4. Vytvorte tabuľku telies Slnečnej sústavy (všetky planéty plus Ceres) s odpovedajúcimi vzdialenosťami (0.38 AU, 0.72 AU, 1.00 AU, 1.52 AU, 2.77 AU, 5.20 AU, 9.54 AU, 19.22 AU a 30.06 AU) a predpovedanými vzdialenosťami podľa Titus-Bodeho rady. Obe vzdialenosti vyneste do grafu a porovnajte.
5. V študijných materiáloch sa nachádza súbor eclipse.txt, ktorý obsahuje tabuľku mesačných zatmení minulého storočia. Prvý stĺpec reprezentuje rok a druhý stĺpec mesiac zatmenia. V treťom stĺpci je typ zatmenia - Čiastočné, Úplné a Polotieňové. Vykreslite si jednotlivé typy zatmení do grafu závislosti na čase a pokúste sa odhadnúť dĺžku sarosu.
6. Z webovej stránky http://www.sidc.be/silso/datafiles stiahnite textový súbor obsahujúci hodnoty priemerného mesačného Wolfovho čísla od roku 1749 (Total sunspot number —y Monthly mean total sunspot number —y TXT). Pozrite si dokumentáciu a vytvorte graf závislosti slnečnej aktivity na čase. Z grafu určte hlavnú periódu slnečnej aktivity.
7. V databáze VizieR nájdite výsledky modelovania lineárneho koeficientu okrajového stemnenia (J/A+A/529/A75/tableu), ktorý je popísaný vzťahom J = J0(l — it(l — cos p)). Z tohto katalógu vypíšte všetky dáta, ktoré zodpovedajú hviezdam hlavnej postupnosti (log g ^ 3.5), namodelované pomocou metódy LSM a modelom ATLAS a súbor uložte. Z neho potom vyberte všetky merania pri použití filtrov U, B a V a vytvorte graf, ktorý bude ukazovať závislosť koeficientu u na teplote s farebne odlíšenými filtrami.
8. Každý deň o 12. hodine bola zaznamenaná pozícia Mesiaca. Nameraná hodnota rektascenzie a deklinácie bola uložená do súboru moon.dat. V prvom stĺpci sa nachádza dátum pozorovania, v druhom je rektascenzia v hodinách a v treťom deklinácia v stupňoch. Určte závislosť denného uhlového pohybu Mesiaca na čase a pomocou vzťahov u = -, rmin = a(l — e) a druhého Keplerovho zákona vypočítajte
1
excentricitu Mesačnej dráhy. Uhlová vzdialenosť medzi dvoma bodmi na sfére o súradniciach a [«2,^2] je
cos A = sinči sin ^2 + cos ô± cos 62 cos(«i — «2).
Výstupom nech je graf závislosti denného pohybu na čase a hodnota excentricity spolu s uvedením vzorca na jej výpočet.
9. Vytvorte simuláciu 2D potenciálu v systéme dvoch korotujúcich telies obiehajúcich okolo spoločného ťažiska po kruhových dráhach. Pomer hmotností telies je 1:10. Potenciál počítajte len v obežnej rovine. Nezabudnite vziať do úvahy aj odstredivú zložku potenciálu! Výsledkom nech je 3D graf, kde osi x a y budú reprezentovať súradnice bodu v obežnej rovine a na ose z bude hodnota potenciálu.
2
Kapitola 2
Príklady použitia numerických metód
1. Libračný bod Zeme Ll je bod medzi Slnkom a Zemou, kde je výslednica všetkých síl nulová. Gravitačná sila Slnka je kompenzovaná gravitačnou silou Zeme a odstredivou silou obežného pohybu. Teleso, ktoré sa nachádza v tomto bode, sa neriadi klasickým Keplerovým zákonom, ale v dôsledku dlhodobej stability systému obieha rovnakou uhlovou rýchlosťou ako Zem. Zo silovej rovnice určte iteračnou metódou vzdialenosť tohto bodu od Zeme v km, ak Ms = 333 000 Mz, r = 149 597 900 km.
2. Halleyova kométa obieha okolo Slnka po extrémne eliptickej dráhe s excentricitou 0.967. Veľká polos jej dráhy je 17.8 AU. Vypočítajte strednú anomáliu M, excentrickú anomáliu E a pravú anomáliu v. V polárnom grafe znázornite polohy na dráhe v ročných intervaloch po prechode perihéliom. Pre výpočet sú potrebné rovnice:
3. Vytvorte funkciu, ktorá bude počítať dynamickú paralaxu dvojhviezd. Vstupné údaje nech sú m a, m b, a a P, ako zdanlivé hviezdne veľkosti oboch hviezd, veľká polos v oblúkových sekundách a obežná doba v rokoch. Funkcia zo vstupných hodnôt s predpokladanými hmotnosťami oboch hviezd rovnými 9JÍ0 a použitím tretieho Keplerovho zákona určí veľkú polos a z nej predpokladanú hodnotu paralaxy. Z modulu vzdialenosti sa vypočítajú absolútne hviezdne veľkosti oboch zložiek, z nich sa určí žiarivý výkon a následne sa vypočítajú približné hmotnosti zo vzťahov hmot-nosť-svietivosť. Tie sa vložia ako vstupné hodnoty do ďalšieho cyklu a celý postup sa opakuje, kým sa nedospeje ku konštantným hodnotám. Výstupom funkcie budú hmotnosti oboch zložiek a hodnota dynamickej paralaxy. Bolometrická hviezdna
2nt
M
M
P ' E — e sin E
r
1 + e cos v
3
veľkosť Slnka je 4.75 mag.
ď_ _ G{mA + gftg)
P2 ~ 4^2
M = m — 5 log r + 5
— = 0.23   — SOT < O.4393t0
é=i-5(S* 2^<^<2°^
-p- = 3200-^-      SDT > 2O93t0
4. Podľa postupu zo skript Obecné astronomie, str. 105-107 vytvorte graf analemy, t.j. závislosť rozdielu pravého a stredného slnečného času na deklinácii Slnka.
4
Kapitola 3
Použitie integračných a diferenciálnych metód
1. Pokúste sa o simuláciu svetelnej krivky tranzitu exoplanéty. Nech má exoplanéta veľkosť 5-násobne menšiu ako materská planéta a nachádza sa od nej v dostatočnej vzdialenosti, aby sme dráhu počas tranzitu pokladali za priamku. Vytvorte simuláciu pre prípad, že planéta prechádza stredom disku a lineárny koeficient okrajového stemnenia je rovný 0.5, pričom platí J = J0(l — it(l — cos p)), kde u je požadovaný koeficient a p je uhlová vzdialenosť od centra disku.
2. Na popis gravitačného potenciálu diskových galaxií sa používa ako jeden z modelov tzv. Toomre-Kuzminov model. Ten využíva axiálnej symetrie pričom predpokladáme, že hmota je rozložená v zanedbateľné tenkom disku, ktorého plošná hustota je definovaná nasledovne:
E (r)
M 2na2
np \ 2
a
-3/2
■y \ 2
a
i -3/2
Vytvorte simuláciu gravitačného potenciálu galaxie v rovine yz. Nech sú všetky konštanty potrebné pri výpočte jednotkové a disk galaxie má polomer R = 2.
3. Železný meteoroid tvaru gule o polomere 10 cm padá zvislým smerom k Zemi. Vo vzdialenosti 500 km od povrchu má rýchlosť 1000 m/s. Vytvorte grafy závislosti polohy, rýchlosti a zrýchlenia na čase, pričom porovnáte tieto veličiny s prípadom bez započítania odporu atmosféry. Odpor vzduchu je popísaný rovnicou Fq = \CSpv2, pričom hustota vzduchu sa mení s výškou a je približne daná vzťahom p = poexp(—-jj), kde po = 1.24 kgm~  a H = 9150 m. Pri výpočte budú potrebné
nasledovné konštanty: Mz = 5.974-1024 kg, G = 6.674- ÍO"11 m3 kg R = 6378 km, pFe = 7870 kgm~3. Abláciu neuvažujte.
C = 0.47,
x = x0 + v0 dt + —a dt v = vq + a dt
5
Kapitola 4
N*"
Statistika
1. a) Zo studijných materiálov načítajte súbor m44.dat, ktorý obsahuje časť výstupu z
fotoaparátu, uložte ho do matice a vykreslite si rozloženie intenzity jednotlivých pixelov pomocou 3D grafu.
b) Vypočítajte aritmetický priemer, medián, modus, štandardnú odchýlku, šikmosť
a špicatosť dátového súboru. Použite definície šikmosti a špičatosti zo skript.
c) Z dátového súboru odstráňte všetky hodnoty vzdialené od mediánu o viac ako
3.5o> [ar = 1.482 median(|x — x\)) pravdepodobne zodpovedajúce hviezdam a znovu zopakujte výpočty z časti b).
d) Určte hodnotu intenzity pozadia a hodnotu šumu.
e) Vytvorte funkciu na výpočet mediánu.
f) Vytvorte skript/funkciu na výpočet modusu.
g) Z pôvodnej matice vytvorte nový dátový súbor, kde budú oblasti šumu nahra-
dené hodnotou pozadia a tieto údaje znovu vykreslite do grafu.
h) Vytvorte nástroj na výpočet inštrumentálnej hviezdnej veľkosti (spolu s veľkosťou
chyby) pre hviezdu s centrom o súradniciach [x, y] a polomere r.
i) Okalibrujte inštrumentálnu hviezdnu veľkosť, ak viete, že 7 Cnc má jasnosť
4.65m. Určte zdanlivú hviezdnu veľkosť hviezdokopy M44. Zistite, aká je limitná hviezdna veľkosť hviezd na snímku (minimálna hodnota S/N = 5).
2. V študijných materiáloch sa nachádzajú dva dátové súbory hdl88041.dat a hr7575.dat. V ôsmom stĺpci matice sa nachádzajú namerané hodnoty hviezdnych veľkostí v rôznych časových okamihoch. Zadefinujte si vlastným spôsobom vytvorenie histogramu (bez použitia vstavanej funkcie numpy.histogram) a aplikujte ho na rozdelenia jasností pre obe hviezdy. Ktorá z týchto dvoch hviezd je premenná a prečo? Svoje rozhodnutie podporte výpočtom indexov normality:
3. Hviezdy A, B a C majú každá inú svetelnú krivku. Hviezda A má konštantnú jasnosť rovnú nule, jasnosť hviezd B a C sa mení v závislosti na fáze <fi nasledovne:
m(B)
m (C) = exp
6
Vytvorte simuláciu pozorovaní týchto hviezd tak, že fázu rovnomerne pokryjete napr. 1000 bodmi a každému bodu vygenerujete hodnotu podľa danej funkcie. Táto hodnota nech je však zaťažená simulovanou nepresnosťou merania, ktorá zodpovedá normálnemu rozdeleniu o pološírke 0.1. Simulované merania vyneste do grafu znázorňujúceho svetelnú krivku a ku každému grafu vytvorte histogram jasnosti. Pre každú situáciu určte hodnoty šikmosti a špičatosti a výsledky diskutujte.
7
Kapitola 5 Lineárna regresia
1. a) Predpokladajte, že všetky hviezdy majú rovnakú jasnosť a ich priestorové rozloženie
je homogénne. Určte koľkokrát sa zvýši pozorovaný počet hviezd zvýšením pozorovacej schopnosti o jednu magnitúdu.
b) Súbor list-ofstars.dat obsahuje zoznam najjasnejších hviezd do 3.0 mag. Na základe odvodenej závislosti určte vhodnú modelovú funkciu, ktorá popisuje vzťah medzi hviezdnou veľkosťou a poradím. Pomocou metódy najmenších štvorcov preložte dáta a určte parametre fit u. Diskutujte rozdiel medzi hodnotou fit u a teoretickou hodnotou z časti a).
2. a) Z katalógu premenných hviezd ASAS vyberte vhodné dáta pre polodotykovú
zákrytovú dvojhviezdu V2509 Sgr. Pomocou uvedenej periódy vypočítajte príslušné fázy a merania hviezdnej veľkosti MAG_3 vykreslite do fázovej svetelnej krivky.
b) Dáta preložte harmonickým polynómom vhodného stupňa. Na výpočet váh
použite chyby merania. Môžte takisto vykresliť samotné chyby príkazom er-rorbar.
c) Predpovedajte hodnoty veličín x2, x\ a s- Následne tieto hodnoty určte výpočtom.
d) Vypočítajte hodnoty ôyp a ôb a vykreslite ohraničenie presnosti fit u yp ± ôyp.
e) Odstráňte odľahlé body a predchádzajúci postup niekoľkokrát zopakujte.
3. V katalógu Vízíer vyhľadajte súbor cefeíd z veľkého Magellanovho mračna B VRI observatíons of LMC Cepheíds, Sebo (2002). Vykreslite závislosť medzi periódou a hviezdnou veľkosťou vo farbe V. Preložte závislosť logaritmickou funkciou s použitím odhadovaných chýb a preveďte ju na vzťah absolútna hviezdna veľkosť-perióda ak viete, že vzdialenosť LMC je 50 kpc.
4. V katalógu Vizier vyhľadajte súbor meraní radiálnych rýchlostí hviezdy HD110555 (Radial velocítíes of HD 110555 (Torres+, 2007)). Vykreslite závislosť radiálnej rýchlosti HRVAa a HRVAb na fáze. Závislosti preložte vhodnou modelovou funkciou. Pri výpočte použite chyby vo forme váh. Z informácií obsiahnutých v dokumentácii a v samotných meraniach skúste odpovedať na otázky:
a) Aký je dvôvod existencie dvoch rôznych fáz PhasA a PhasB? Ako funguje tento
hviezdny systém?
b) Aká je približne hodnota vzdialenosti zložiek A, ak hmotnejšia z nich má pri-
bližne hmotnosť 1 M0?
8
Kapitola 6 Nelineárna regresia
1. Súbor cont0109a.txt obsahuje spektrum v oblasti absorbčnej vodíkovej čiary H^. Na fitovanie spektrálnej čiary použijeme Gaussovu a Lorentzovu funkciu:
(x-xa)2
--J-
G = G0 + Ae 2CT L = L0 + A
1 (x-x0)2+j2
Vykreslite si spektrum do grafu, zvoľte počiatočné hodnoty parametrov a preložte dáta oboma krivkami. Pomocou metódy x2 vyhodnoťte, ktorá funkcia je vhodnejšia. Akú vlnovú dĺžku má centrum čiary?
2. Z databázy mCPod vyberte všetky pozorovania hviezdy HD 81009 v Strômgrenovom filtri v. Prvý stĺpec zodpovedá označeniu hviezdy, druhým je HJD, v treťom stĺpci je jasnosť a v šiestom indentifikačné číslo pozorovateľa. Perióda hviezdy je 33.978462 dňa. Ako modelovú funkciu použite
/ = Í2Ó2 + + f37ô37 + /noíno + A exp(l - COsh ^0-round(^o)))
kde ôq je rovné 1 v prípade, že meranie vykonal pozorovateľ s identifikáciou q a 0 inak. Parametre funkcie sú f2, f27, f37, f 110, A, Lpo a a. Na začiatku si vytvorte graf, z ktorého určíte počiatočné hodnoty parametrov. Vykonajte nelineárnu regresiu s jednotkovými váhami. Výsledkom nech je graf s dátami posunutými na nulovú hladinu (závislosť y —      fqôq na fáze).
3. Na webových stránkach ESO sa nachádzajú spektrá niekoľkých desiatok hviezd. Stiahnite si spektrum hviezdy HD22049 spektrálneho typu K2V a pomocou Planc-kovho zákona určite jej teplotu. Doporučený postup:
a) Vykreslite si spektrum hviezdy. Definujte kontinuum a to tak, že si súbor rov-
nomerne rozdelíte na intervaly vlnových dĺžok (napr. po 10 nm) a z každého intervalu vyberiete meranie s maximálnou hodnotou žiarivého toku.
b) Teoretické rozloženie toku je dané Planckovým vzťahom: F = ^exp^b/x)-i^ ^e
B = jr^f- Odhadnite približnú hodnotu parametra B a dopočítajte odhad parametra A. Ako vstupné hodnoty použite len vybrané hodnoty z časti a), na ktoré aplikujete nelineárnu regresiu.
c) Z hodnoty parametra B vypočítajte približnú povrchovú teplotu hviezdy.
9
Kapitola 7
Periodogramy, O-C diagramy
1. a) Conroy et al. (2014) publikoval katalóg presných pozícií zákrytov tranzitujúcich
premenných hviezd na základe meraní družice Kepler (Kepler mission. IV. Eclípse times f or close bínaríes (Conroy+, 2014))- Niektoré systémy však vykazujú systematické odchýlky od vypočítaných hodnôt času predpokladaných zákrytov. V databáze ViZieR vyhľadajte príslušný katalóg a vyberte z neho dáta ku hviezde KIC11247386.
b) Vykreslite si časovú zmenu rezíduí ETV, odhadnite hodnotu periódy a polohu
MO a následne vypočítajte príslušné fázy.
c) Závislosť ETV na fáze preložte funkciou f = fo + Acos[2tt(<f) — 0o)]- Počítajte s
váhami nepriamo úmernými druhej mocnine chýb.
d) Postup opakujte pre rovnomerne pokrytý interval vhodných periód. Vytvorte
závislosť štandardnej odchýlky na perióde (periodogram) a určte skutočnú hodnotu periódy.
2. V študijných materiáloch sa nachádza dátový súbor neznámej hviezdy. V prvom stĺpci je heliocentrické juliánske dátum, v druhom stĺpci hviezdna veľkosť v mag-nitúdach. Pomocou lineárnej regresie za použitia harmonických polynómov druhého rádu vytvorte periodogram tejto hviezdy, ak predpokladáte, že svetelné zmeny sú dané rotáciou hviezdy hlavnej postupnosti. Pokúste sa určiť správnu rotačnú periódu neznámej hviezdy.
3. a) SQ Tau je zákrytová dvojhviezda, ktorej zložky sú horúca hviezda triedy B
a červený M trpaslík. Keďže pomer jasností týchto zložiek je veľmi veľký, svetelná krivka obsahuje len primárne minimum. Perióda svetelných zmien je 230.01323 dňa a určená hodnota nulového minima MO = 42446.22 d. Dátový súbor SQ.dat obsahuje v prvom stĺpci hodnotu MJD a v druhom jasnosť vo filtri V. Vykreslite si závislosť jasnosti na MJD.
b) Preložte závislosť modelovou funkciou / = /o + X^ ^exP (~ ^~2°ŕ ) • Počiatočné
hodnoty t0i určte z polohy jednotlivých miním. Váhy neuvažujte.
c) Z hodnôt P a MO určte „calculated" hodnoty polôh miním, ztotožníte ich s
„observed" hodnotami t0i a určte rozdiel O-C. Vypočítajte príslušné epochy a vykreslite O-C diagram, kde budú znázornené aj chybové úsečky hodnôt miním (ôtoi). Vzniknutú závislosť preložte vhodnou funkciou a správanie periódy správne interpretujte.
10
4. V databáze 0-C brána vyhľadajte dáta ku hviezde RT And a uložte ich do dátového súboru. Z neho potom vhodným spôsobom načítajte prvé tri stĺpce bez toho, aby ste pôvodný súbor museli upravovať. Vykreslite si O-C diagram, preložte ho vhodnou modelovou funkciou a interpretujte charakter periódy.
11
Kapitola 8 Bootstrap technika
1. a) Dátový súbor data.txt obsahuje pozorovania okamihu tranzitu exoplanéty. V
hviezdna veľkosť. Z modelovania svetelnej krivky takéhoto zákrytu je možné určiť rôzne parametre, ako napríklad polomery exoplanéty a materskej hviezdy. Okrem hodnôt samotných sú dôležité aj veľkosti ich neistôt. Tie sú však pri použití nelineárnej regresie často skreslené. Jeden zo zjednodušených mode-
b) Použite techniku bootstrap, definovanú nasledovným postupom. Vytvorte množinu nových pozorovaní tak, že z pôvodných dát zakaždým náhodne vyberiete n pozorovaní s možnosťou opakovania, pričom n je zhodné s počtom pôvodných pozorovaní. Na novej množine zakaždým zopakujte proces nelineárnej regresie a uložte si výsledky parametrov. Tento cyklus zopakujte 1000-krát, čím získate sadu 1000 vektorov parametrov. Vypočítajte hodnoty štandarných odchýlok jednotlivých parametrov a porovnajte ich s hodnotami neistôt z časti a).
prvom stĺpci sa nachádza modifikovaný juliánsky dátum, v druhom zdanlivá
12