1
[0] ÚVOD: [0.1]
Stručná historie
Zde bych rád napsal několik málo řádek o budování systému GPS NAVSTAR. To začalo v roce 1973 a bylo koncipováno jako obranný navigační sytém Spojených Států Amerických. Vedením tohoto
projektu
bylo
pověřeno
U.S.
Air
Force
(Letectvo
Spojených Států). Dále s touto institucí spolupracovalo U.S. Army
Navy
(Vojenské
námořnictvo
Spojených
Států)
a
DMA
(Defense Mapping Agency). O pět let později, v roce 1978, se k tomuto projektu připojilo i dalších devět členských států NATO (North Atlantic Treaty Organization). Požadavky na vybudování takového navigačního systému byly, aby
bylo
možno
v reálném
čase
kdykoliv
a
kdekoliv
zjistit
přesné určení polohy (řádově do 10 m) libovolného počtu i rychle se pohybujících objektů.
[0.2]
Systém GPS-NAVSTAR
Systém GPS-NAVSTAR (Global Positioning System – NAVigation System
using
Time
And
Ranging)
tvoří
tři
části:
řídící,
kosmická a uživatelská. Řídící část tvoří sledovací stanice (obr. č.:1) rozmístěné po celé Zemi a hlavní řídící stanice (Colorado Springs), jež
obr. č.: 1
2
zpracovává telemetrické informace a výsledky sledování pohybu družic z ostatních sledovacích stanic a přes jejich vysílače tyto informace předává jednotlivým družicím spolu s povely pro řízení provozního režimu a korekcemi drah. Kosmickou nominálně
část
24
družic
rozmístěných
na
tvoří rovnoměrně
šesti
oběžných
drahách (obr. č.:2). Dráhy družic mají
sklon
kruhové. povrchem
55o
a
Výška Země
je
jsou
téměř
družic
nad
20
183
km
a
jejich oběžná doba činí 11 hod 58 min,
tedy
12
hvězdných
hodin.
Družice (obr. č.:3) mají hmotnost 845 kg a přibližné rozměry jsou
obr. č.: 2
2,0 x 1,0 x 1,5 metrů. Životnost družic je počítána na 7,5 roku provozu. Palubní baterie jsou dobíjeny
slunečními
původně
zabudovávaly
články
o
ploše
rubidiové
m 2.
7,15
oscilátory,
Do
které
družic jsou
se
dnes
nahrazovány oscilátory vodíkovými či cesiovými.
Tyto
oscilátory
mají
za
úkol udržovat velmi přesné časové a kmitočtové informace. A
nakonec
uživatelský
segment.
Ten je tvořen přijímači signálů GPS pomocí obr. č.: 3
zařízení.
antén Signál
koherentních
a
registračních
GPS
tvoří
řada
kmitočtů,
které
jsou
odvozeny ze základní frekvence f0 = 10,23 MHz. Ta je ale kvůli kompenzaci průměrného relativistického efektu snížena o 4,45 x 10-10
f0. Tento základní kmitočet je udržován oscilátory na
družicích s relativní přesností lepší než 10-13.
3
[0.3]
Použitá symbolika v kapitolách [1.5] - [1.7]
p, p
vektor, skalár (velikost vektoru)
•
skalární součin
⊗
vektorový součin
Fj
Delauneyovy proměnné, j=(1, 2, …, 5)
z, θ, ζ
precesní úhly
Δψ , Δε
nutace v délce, nutace ve sklonu
ε0
střední sklon ekliptiky
g
střední anomálie Země
L , L1, L2, L3
vektorový integrál ploch, jeho složky
H
integrál energie, hamiltonián
λ , λ 1, λ 2, λ 3
vektorový integrál Laplaceův, jeho složky
T, V
kinetická energie družice, její silová funkce
GM
geocentrická gravitační konstanta
[ρ,u]
polární souřadnice družice v rovině její dráhy
m
hmotnost družice
ρ , x1 , x2 , x3
geocentrický průvodič družice, jeho složky
r , x1 , x2 , x3
geocentrický průvodič družice, jeho složky
x , x1 , x2 , x3
vektor polohy družice, jeho složky
ρ=r=x
platná relace
v , v1, v2, v3
vektor rychlosti družice, jeho složky
I
výstupní uzel dráhy družice
Ω
rektascenze výstupního uzlu dráhy družice
e
excentricita dráhy družice
i
sklon dráhy družice
a
velká poloosa dráhy družice
ν
pravá anomálie dráhy družice
E
excentrická anomálie družice
ω
argument perigea dráhy družice
u
argument deklinace dráhy družice
4
[1] TEORETICKÉ ZÁKLADY: [1.1]
Fourierova transformace
Pro interpolaci dat, kterou jsem v této práci prováděl, bylo nutné nejprve určit frekvence (resp. periody) a amplitudy harmonických
průběhů
terestrickými elementy. periody)
hodnot
souřadnicemi
V některých lze
transformace
vyčíst
ať
nebo
přímo
případech
tyto
přímo
představuje
už
Keplerovými hodnoty
z vynesených
mnohem
reprezentovaných (amplitudy
dat,
přesnější
dráhovými
nicméně
nástroj
jak
a
tato se
těchto hodnot dopátrat. Analýzu výchozích dat jsem prováděl v systému Matlab 4.0, který tuto transformaci po výpočetní stránce zajišťoval. Zde bych
chtěl
jenom
nastínit
Fourierovu
analýzu
po
teoretické
stránce. Mějme tedy definovanou funkci g(x), která splňuje podmínku absolutní integrovatelnosti: ∞
∫ g ( x) dx < ∞
−∞
a
vyhovuje
i
ostatním
Dirichletovým
podmínkám.
Tuto
funkci
chceme rozvinout do řady Fourierových funkcí: 1, cos(x), sin(x), cos(2x), sin(2x), … což je systém funkcí, který je na intervalu <-π,π> ortogonální a tvoří bázi. Funkci g(x) tedy můžeme rozepsat do Fourierovy řady:
g ( x) =
∞ a0 + ∑ [a k . cos(kx) + bk . sin(kx)] 2 k =1
5
přičemž koeficienty ak a bk mají tvar:
ak = bk =
1
π 1
π
π
∫ g ( x). cos(kx).dx
−π
π
∫π g ( x). sin(kx).dx
−
a říká se jim Fourierovy koeficienty. Ty vypočteme na základě MNČ minimalizací výrazu: π
g ( x) − g N ( x) =
∫π (g ( x) − g
N
)
( x) 2 dx
−
v němž funkce gN(x) má tvar:
g N ( x) =
N a0 + ∑ [a k . cos(kx) + bk . sin(kx)] 2 k =1
Za použití Eulerových vzorců pro přechod ke komplexním číslům se rozvoj funkce g(x) ještě zjednoduší. Tyto vzorce mají tvar: e ix + e − ix 2 ix e − e −ix sin( x) = 2i cos( x) =
a původní báze se změní na bázi komplexních čísel: 1, eix, e-ix, e2ix, e-2ix, … Potom funkci g(x) můžeme psát ve tvaru: g ( x) =
1 2π
∞
∑c
k = −∞
k
.e ikx
Koeficienty ck mají pak tvar:
ck =
1 (a k − i. sgn(k ).bk ) nebo jinak: c k = 1 2π 2
π
∫π g ( x).e
− ikx
.dx
−
Až do této chvíle jsme se pohybovali pouze v intervalu <-π,π>. Nyní tento interval rozšíříme.
A to
substitucí
τ∈<-N,N>. A dále:
τ=x
N
π
a po derivaci dx =
N
π
dτ
proměnou
6
Změní se také Fourierovy koeficienty, které nabudou tvaru: N
1 ck = 2N
∫ g (τ ).e
− ik
π N
τ
.dτ
−N
a rozvoj funkce g(x) bude vypadat takto:
g ( x) =
Logický
požadavek
je,
⎡ 1 ⎢ ∑ k = −∞⎣ 2 N ∞
N
∫ g (τ ).e
− ik
π N
τ
−N
aby
se
N
⎤ ik π x .dτ ⎥.e N ⎦
Æ ∞. Proto zavedeme tuto
symboliku:
1 k = k .Δf = f k = Δf a dále 2N 2N Tím rozvoj funkce g(x) přejde na tvar:
⎡N ⎤ g ( x) = ∑ ⎢ ∫ g (τ ).e − 2πif kτ dτ ⎥.e 2πif k x .Δf k = −∞⎣ − N ⎦ ∞
a pakliže vyjdeme z našeho požadavku aby N Æ ∞, musí Δf Æ 0 a dále sumační znak přejde na integrál a fk přejde na f. Tedy:
⎡∞ ⎤ 2πifx − 2πifτ ∫−∞⎢⎣−∫∞g (τ ).e .dτ ⎥⎦.e .df ∞
g ( x) =
Výraz
v hranatých
závorkách
označíme
jako
G(f).
Takto
dostáváme dvě funkce: ∞
g ( x) = ∫ G ( f ).e 2πifx .df −∞
∞
G( f ) =
∫ g ( x).e
− 2πifx
.dt
−∞
Funkce dat.
G(f)
Funkce
transformace.
nazýváme g(x)
je
Originál
Fourierovou dvojicí.
Fourierovým pak
obrazem
nazývána
(g(x))
a
jako
obraz
nebo
též
zpětná
(G(f))
spektrem
Fourierova
nazýváme
též
7
Problémem Fourierovy transformace je ale to, že tato není schopna detekovat takové periodické jevy, které mají frekvenci blížící
se
počtu
transformace,
měření.
ve
smyslu
Nejlepších
výsledků
důvěryhodnosti,
se
této
dosahuje
v případech, kdy frekvence periodických jevů je menší než 1/2 počtu měření.
[1.2]
Harmonická anlýza Fourierovou řadou
Používá
se
vykazují-li
data
jasnou
periodicitu.
Potom
považujme interval od bodu x=a do bodu x=b, za nímž se průběh naměřených
dat
alespoň
zhruba
opakuje,
za
periodu
P=b-a.
Fourierova řada bude mít tedy tvar:
yi = A0 + a.sin (t i + A) + b.sin (2t i + B ) + c.sin (3t i + C ) + ... přičemž
t
je
nová
proměnná,
která
se
získá
z
proměnné
x
transformací:
ti =
2π xi P
V obou posledně uváděných vztazích je index i={1,2,3,…,n}, kde n značí počet dat vstupujících do vyrovnání a je tedy shodný s počtem
rovnic,
ze
kterých
jsou
pak
následně
počítány
koeficienty interpolační funkce. Je dále výhodné Fourierovu řadu upravit (zjednoduší se derivace Fourierovy řady) pomocí substituce:
a.sin A = A1 a. cos A = A2
,
b. sin B = B1 b. cos B = B2
,
c. sin C = C1 c. cos C = C2
, …
na tvar:
yi = A0 + A1 cos t i + A2 sin t i + B1 cos 2t i + B2 sin 2t i + C1 cos 3t i + C 2 sin 3t i + ... Rovnice oprav tedy bude vypadat takto:
vi = A0 + A1 cos t i + A2 sin t i + B1 cos 2t i + B2 sin 2t i + C1 cos 3ti + C 2 sin 3t i + ... − yi
8
Úloha tak přejde na vyrovnání zprostředkujících měření, kde neznámé parametry jsou A0, A1, A2, B1, B2, C1, C2, …. Členy matice plánů A (design matrix A) pak tvoří parciální derivace Fourierovy harmonické řady:
mat Ai , j = def{ 1, cos t i , sin t i , cos 2t i , sin 2t i , … }i , j Vektor neznámých je tedy dán vztahem:
h = h 0 + dh kde h0 je vektor přibližných neznámých, vektor dh je roven výrazu:
(
)
−1
dh = AT . A . AT .l ' , který je dán podmínkou:
vT.v Æ min
přičemž vektor l’ je vektor redukovaných měření a je roven:
l ' = l − l 0 a l 0 = A.h 0 a to za předpokladu, že vektor oprav jednotlivých měření (v našem případě souřadnic resp. Keplerových dráhových elementů) je: v = A.h − l
A dále zavedeme kovarianční matici neznámých parametrů:
Q = (AT . A)
−1
potom střední chyby neznámých parametrů Fourierovy řady jsou dány vztahem:
mhi = m0 . Qi ,i Výraz m0 je odhad jednotkové střední chyby:
m0 =
v T .v n−k
9
kde n je počet všech měření (v podobě jednotlivých souřadnic nebo
Keplerových
dráhových
elementů)
a
k je
počet
nutných
měření.
[1.3]
ITRS a rámec ITRF
Souřadnice družic GPS tak jak jsou uváděny a publikovány v oficiálních
datových
souborech
jsou
v systému
ITRS
(International Terrestrial Reference System). Tento systém je realizován souřadnicemi bodů ITRF (International Terrestrial Reference Frame). Počátek tohoto referenčního rámce (frame) je definován
v těžišti
Země
a
to
s přesností
určení
10
cm.
Referenční pól (je identický s osou Z kartézské soustavy) a referenční poledník (je identický s rovinou tvořenou osami XZ) jsou
konzistentní
Terrestrial
s odpovídajícími
System,
s přesností
BIH
0,003“.
–
Bureau
Osa
X
směry
systému
International leží
BTS
de
v nultém
(BIH
l’Heure) poledníku
procházejícím Greenwichem a v rovině rovníku. BIH referenční pól
byl
připojen
k CIO
(Conventional
International
Origin)
s přesností 0,03“. Vzhledem k pohybům tektonických desek, jsou časové změny souřadnic
bodů
definujících
rámec
ITRF
určovány
kombinací
přímého měření a teoretických hodnot vyšlých z modelů pohybů a to za podmínky tzv. NNR (No Net Rotation). Tj. za podmínky nulové
rotace
minimálního
sítě,
což
odpovídá
diferenciálního
Tisserandově
přírůstku
podmínce
kvadrátu
momentu
hybnosti. Protože tento systém je pevně spjat s rotující Zemí, není systémem inerciálním.
[1.4]
ICRS a rámec ICRF
Vzhledem
k tomu,
že
Newtonův
gravitační
zákon
platí
v inerciálním systému a pohybové rovnice jsou v něm jednodušší (nemusí se zavádět korekce z neinerciality systému – zrychlení Eulerovo,
Coriolisovo,
dostředivé)
je
výhodné
pro
výpočet
10
Keplerových
dráhových
elementů
vyjít
právě
z onoho
inerciálního systému. Takovým systémem je ICRS (International Celestial souborem tzv.
Reference
System).
mimogalaktických
kvasarů,
které
Systém
zdrojů tvoří
ICRS
je
realizován
elektromagnetického referenční
rámec
záření, –
ICRF
(International Celestial Reference Frame) . Počátek ICRF je v barycentru sluneční soustavy. Osa Z kartézského souřadného systému
směřuje
do
polohy
středního
rotačního
pólu
v epoše
J2000.0. Osa X směřuje do Jarního bodu v epoše J2000.0 a osa Y doplňuje
tento
systém
na
pravoúhlý
v matematicky
kladném
směru. Zde bych chtěl uvést na pravou míru to, proč jsem jako inerciální systém označil systém ICRS. Pravda je ta, že přísně inerciální systém být ICRS nemůže a to proto, že je vázán na vesmírné objekty, které vykonávají svůj vlastní pohyb. Nicméně tento pohyb je vzhledem k obrovským vzdálenostem kvasarů tak malý,
že,
přihlédneme-li
k naší
přesnosti
úhlového
měření,
chyby způsobené tímto zjednodušením jsou zanedbatelné.
[1.5]
Transformace [x,y,z]ITRS Æ [x,y,z]ICRS
Abychom mohli spočítat Keplerovy dráhové elementy družic je nejprve nutné provést transformaci terestrických souřadnic družice
na
její
hvězdné
souřadnice.
Tato
transformace
je
realizována maticí rotace, které sestává z pěti dílčích matic rotace. Ty popisují vliv precese, nutace, rotace Země a pohybu pólu na terestrické souřadnice. Samotná transformace je dána vztahem: R = RPT .RNT .RTT .R Xp .RYp
a tedy ⎛X⎞ ⎛X⎞ ⎜ ⎟ ⎜ ⎟ = R.⎜ Y ⎟ ⎜Y ⎟ ⎜Z⎟ ⎜Z⎟ ⎝ ⎠ ICRS ⎝ ⎠ ITRS
11
[1.5.1]
Vliv precese
Precese je pohyb rotační osy Země okolo pólu ekliptiky a dále pohyb tohoto pólu vůči inerciální soustavě realizované rámcem ICRF. Je to pohyb téměř věkovitý (sekulární) a dá se rozdělit na složku
planetární
způsobena složka,
Měsícem jak
způsobena
již
a
lunisolární.
a
Sluncem. název
působením
Lunisolární
složku
napovídá,
planet. nazveme
je
Planetární je
Složku
lunisolární nazveme lunisolární precesí a planetární
složka
analogicky
ε = 23,5o Pe
planetární precesí. Působením obrázek
č.:4)
lunisolární se
střední
precese
(viz.
světový
pól
PN obr. č.: 4
pohybuje okolo pólu ekliptiky. Tento pohyb má periodu 25 800 let a říká se jí Platónský rok. Poloměr tohoto pohybu je roven sklonu ekliptiky vůči rovníku, a to 23,5o. Proto se mění i poloha jarního bodu a to o 50,3“ za rok. Planetární složka precese má pak za následek pohyb pólu ekliptiky vůči inerciální soustavě. Pohyb pólu ekliptiky je téměř
kruhový s poloměrem 90’ a periodou
70000let. Vliv
planetární precese je mnohem menší než lunisolární. Precesní matice vyjadřující působení precese je součinem tří dílčích rotačních matic. Pro naplnění dílčích rotačních matic je třeba spočítat hodnoty precesních úhlů. Tyto úhly se počítají řadami, které jsou níže uvedeny. V těchto řadách se počítá s časovým intervalem T mezi danou epochou a epochou J2000.0
v Juliánských
stoletích.
Při
výpočtu
T
by
bylo
správnější použít barycentrický dynamický čas (TBD). Ale pro naše
účely
postačí
počítat
s terestrickým
dynamickým
časem
(TDT). Ten je s barycentrickým časem, až na periodické variace vůči času atomovému, shodný. Jeho výpočet je také jednodušší.
12
Výpočet terestrického dynamického času je dán vztahem:
TDT = TGPS +
19 + 32,184 86400
TGPS se udává v MJD (Modifikované Juliánské Datum) a pak i TDT vyjde v MJD. Výpočet precesních úhlů je dán řadami:
z=
2306,2181. T + 1,09468. T 2 + 0,018203. T 3 206264,8
θ=
2004,3109. T − 0,42665. T 2 − 0,041833. T 3 206264,8
ς=
2306,2181. T + 0,30188. T 2 + 0,017998. T 3 206264,8
kde T je časový interval mezi danou epochou a epochou J2000.0 vyjádřen v juliánských staletích o délce 36525 dnů a je tedy dán vztahem:
T=
TDTMJD − 51544,5 36525
Precesní matice má tvar:
R = Rz (− z ).R y (θ ).Rz (− ζ )
[1.5.2]
Vliv nutace
Nutační pohyb je pohyb, který vykonává pravý (okamžitý) světový pól okolo středního světového pólu. Tento pohyb je způsoben
periodickými
změnami
polohy
Slunce
a
Měsíce
vůči
Zemi. Hlavni nutační perioda je 18,62 roku a její amplituda má hodnotu 9,21”. Schematicky vypadá nutační pohyb jak je ukázáno na obrázku č.:5.
13
Nutace
má
dvě
složky,
ε = 23,5
rozdíl
o
15,6’
PN 9,21“
Δψ
v délce
a
nutaci
ve sklonu Δε. Nutace v délce je
dráha středního světového pólu
Pe
nutaci
mezi
ekliptikální
délkou
ovlivněnou
délkou
nutací ve
nutací
a
neovlivněnou.
dráha pravého světového pólu
Nutace
sklonu
PN
vůči rovině rovníku ovlivněné
rozdíl
mezi sklonem roviny ekliptiky nutací
obr. č.: 5
a
rovině
neovlivněné
nutací.
rozlišujeme
ještě
dlouhoperiodickou a Δε’) a krátkoperiodickou
je
rovníku Dále nutaci
(složky
Δψ’
(složky dψ a dε).
Nutace v délce a nutace ve sklonu rozepsaná podle tohoto rozdělení mají tvar:
Δψ = Δψ’ + dψ Sklon
roviny
ekliptiky
vůči
a
Δε = Δε’ + dε
rovině
rovníku
(střední
sklon
ekliptiky) je dán řadou:
ε 0 = 23°26'21,448"−46,8150". T − 0,0059". T 2 + 0,001813". T 3 kde je T časový interval mezi danou epochou a epochou J2000.0 vyjádřen v juliánských staletích o délce 36525 dnů a je tedy definován vztahem:
T=
TDTMJD − 51544,5 36525
Další nutační parametry jsou dány řadami: 263 ⎛ 5 5 ⎞ Δψ = ∑ ⎜⎜ ( Ai + A'i .T ).sin(∑ N j F j ) + A"i . cos(∑ N j F j ) ⎟⎟ i =1 ⎝ j =1 j =1 ⎠ 263 ⎛ 5 5 ⎞ Δε = ∑ ⎜⎜ ( Bi + B'i .T ).cos(∑ N j F j ) + B"i .sin(∑ N j F j ) ⎟⎟ i =1 ⎝ j =1 j =1 ⎠
14
kde Fj jsou Delauneyovy proměnné, které postupně znamenají: F1 = l =
střední anomálie Měsíce
F2 = l’ = střední anomálie Slunce F3 = F = L - Ω =
střední délka Měsíce
- střední délka
výstupního uzlu Měsíce F4 = D =
střední elongace Měsíce
F5 = Ω =
střední délka výstupního uzlu Měsíce
Hodnoty Ai, Ai’, Ai”, Bi, Bi’, Bi”, Nj jsou uvedeny v příloze č.: 1. Delauneyovy proměnné jsou dány řadami, které jsou uvedeny v literatuře [2] na straně 101. Nutační matice má tvar:
RN = Rx (− ε 0 − Δε ).Rz (− Δψ ).Rx (ε 0 ) [1.5.3]
Vliv rotace Země – otočení o hvězdný čas
Protože rámec ITRF je připojen
k rotující
Zemi,
je třeba tyto souřadnice od ní
„odpojit“.
li
Přenásobíme-
tyto
příslušnou
souřadnice maticí
rotace,
získáme souřadnice, jejichž osa
Z směřuje
světového
pólu
do
pravého
v příslušné
obr. č.: 6
epoše. Hvězdný čas, o který se bude soustava otáčet je čas GST (Greenwich Sideral Time), tedy Greenwichský hvězdný čas. Čas GST je závislý na časovém rozdílu mezi okamžikem výpočtu a epochou
J2000.0.
Je
ale
nutné
tento
časový
rozdíl
počítat
v čase UT1, protože ten nejlépe odpovídá kolísavosti rychlosti rotace Země, jejíž velikost dosahuje řádově milisekund (viz. obr. č.:6).
15
Předpokladem k tomuto výpočtu je znalost času TGPS. Tedy známe-li TGPS můžeme vypočítat UT1 pomocí tohoto vztahu: UT1 = TGPS–(TGPS–UTC)+(UT1–UTC) Hodnota (UT1–UTC) je udávána v souborech se souřadnicemi pólu
(*.erp).
hodnota
Díky
kolísá.
nepravidelnosti
Proto
je
službou
v rotaci
IERS
Země
i
(International
tato Earth
Rotation Service) udržován tak, aby tato hodnota byla menší než je 1 vteřina. Hodnota (TGPS-UTC) byla do 30.6.1997 včetně rovna: (TGPS-UTC)= 11 sec. Výpočet
Greenwichského
zdánlivého
hvězdného
času
je
dán
vztahem: rad SGA = SGM + Δψ .cos(ε 0 + Δε )
kde výrazy Δψ, Δε, ε0 jsou nutační úhly počítané dle vztahů v kapitole hvězdný
[1.5.2]
čas
a
hodnota
v radiánech.
Ten
rad SGM
se
je
střední
vypočte
dle
Greenwichský následujících
vztahů a řad:
S GM = 24110,54841 + 8640184,812866 . T + 0,093104 . T 2-6,2.10-6 . T 3 + c kde písmenko c představuje korekci času pro UT1 ≠ 0: ⎛ 43200 + 31557610 . 4 .T ⎞ c = 43200 + 315576.104 . T − int⎜ ⎟ .86400 86400 ⎠ ⎝
16
Hodnota
časového
rozdílu
T
(mezi
danou
epochou
a
epochou
J2000.0) je počítána v čase UT1 v juliánských stoletích. Je dána vztahem:
T=
UT 1MJD − 51544,5 36525
A dále je třeba hodnotu SGM převést na radiány:
rad SGM =
SGM .π ⎛ S .π ⎞ − int⎜ GM ⎟ .2π ⎝ 43200.2π ⎠ 43200
Rotační matice otočení o hvězdný čas bude mít tedy tvar:
RT = R z (− S GA ) [1.5.4]
Vliv pohybu pólu – volná nutace
Vektor
pravého
světového
pólu
se
nachází
velmi
blízko
okamžitého vektoru rotace Země. Tato rotační osa mění stále svoji
polohu
v tělese
Země
a
to
způsobuje
stálou
změnu
souřadnic pozorovacích míst na povrchu Země. Euler a později i Chandler stanovili, že pravý světový
pól
přibližně
vykonává
pohyb
kruhový
hlavní
osy
(pohyb
bývá
okolo
setrvačnosti také
volná
nutace).
perioda
tohoto
nazýván Hlavní
pohybu
tzv.
Chandlerova
její
naměřená
je
perioda
hodnota
a
činí
433 dnů. Maximální amplituda pohybu
dosahuje
hodnoty
0,25“. Pohyb rotační osy je znázorněn na obrázku č.:7.
obr. č.: 7
17
Aby nedocházelo ke stálým změnám souřadnic pozorovacích stanovisek na Zemi byl v letech 1900 – 1905 určen průměrný vektor rotace Země a do něj byla vložena osa Z systému ITRS. Takto
stvořený
International
pól
je
Origin).
Od
označován tohoto
jako
vektoru
CIO se
(Convetional dále
počítají
souřadnice pravého (okamžitého) pólu. Tyto souřadnice jsou určovány z měření VLBI (Very Long Baseline Interferometry), SLR (Satellite Laser Ranging) a GPS a bývají uváděny v bulletinech službou IERS. Souřadnou soustavu v níž jsou uváděny souřadnice pravého pólu
tvoří
osa
x
-
je
vložena
do
směru
Greenwichského
poledníku a osa y - je vložena do směru poledníku o západní zeměpisné
délce 90o.
Matice rotace vyjadřující pohyb pólu jsou následující:
R Xp = Rx ( y p ) RYp = R y (x p ) [1.5.5]
Časové systémy použité při výpočtech
Vzhledem
k tomu,
že
jsem
při
řešení
transformace
mezi
terestrickým systémem souřadnic a hvězdným systémem souřadnic použil hned několik časových systémů, považuji za vhodné se o nich zmínit ve zvláštní kapitole. Časy, které jsem v této práci použil, ať již přeneseně nebo přímo, jsou časy rotační, dynamické a atomové. Čas UT1 je časem rotačním a je použit pro výpočet matice otočení o hvězdný čas. A to proto, že, jak již bylo výše zmíněno, tento čas nejlépe vyhovuje nerovnoměrnostem v rotaci Země, způsobeným mnoha fyzikálními vlivy. Sekunda UT1 je kvůli těmto
nerovnoměrnostem
v rotaci
Země
delší
než
sekunda
TAI
(viz. dále), což se projevuje rozbíhavostí času rotačního a času atomového.
18
Čas čas,
TAI
který
v roce
1967
je
atomový
byl
zaveden
a
jehož
jedna
sekunda byla definována na základě záření atomu Cesia Cs133. byl
Počátek určen
shodoval
tak,
času
aby
s časem
1.1.1958. obr.č.: 8
tohoto
se
UT2
Velikost
dne
atomové
sekundy byla stanovena tak,
aby odpovídala efemeridové sekundě v roce 1900. Čas
UT2
greenwichském
je
definován
poledníku
jako
světový
zjištěný
čas
na
astronomickým
okamžitém měřením
a
opravený o pohyb pólu a vlivy sezónních variací v rychlosti rotace Země. Čas UTC, který je použit v rozdílu (UT1-UTC), tabelovaném v souborech souřadnic pólu *.erp, je založen na základě času TAI.
UTC
je
v přibližné (UT1-UTC)
používán
shodě je
v občanském
s časem
udržován
rotačním.
tak,
aby
životě To
a
je
znamená,
nepřekročil
udržován že
rozdíl
hodnotu
1s
přičtením resp. odečtením právě jedné sekundy k resp. od času UTC. Časový rozdíl (UTC-TAI) činí dnes 32 sekund. Dalším časovým systémem použitým při výpočtech je čas GPS, TGPS. Je založen na čase atomovém a byl definován tak, aby se v roce 1980 rovnal času UTC. Jeho rozdíl od času UTC je stálý a má hodnotu 19 sekund. Dále
byl
použit
časový
systém
TDT,
tedy
terestrický
dynamický čas. Ten navazuje na efemeridový čas a platí, že 1.1.1977 0h 0m 0s TAI = 1.1.1997 0h 0m 32,184s TDT. Na základě tohoto vztahu lze tedy spočítat
čas TDT při znalosti
času
atomového. Barycentrický
dynamický
čas
TDB
nebyl
při
výpočtech
použit, protože byl nahrazen časem TDT, který je výpočetně jednodušší a pro naše účely plně postačil. Nicméně i o něm
19
bych
se
rád
zmínil.
Čas
TDB
se
liší
od
času
TDT
pouze
periodickými variacemi, které vyjadřují tyto vztahy:
TDB = TDT + 0,001658 s.sin ( g + 0,0167.sin ( g )) přičemž g je střední anomálie Země a je dána vzorcem:
(
g = 357,528o + 35999,050 o.T
a
T
je
časový
interval
mezi
2π )360
epochou
o
J2000.0
a
počítanou
epochou v juliánských staletích. Obrázek č.:8 přehledně ukazuje vzájemné vztahy časových systémů, které jsem výše popsal.
[1.6]
Výpočet stavových vektorů
Abychom
mohli
z hvězdných
souřadnic
spočítat
Keplerovy
dráhové elementy, je nutné určit nejprve tzv. stavové vektory. Stavový
vektor
definuje
jak
polohu
družice
v daném
časovém
okamžiku, tak i její rychlost. Pro výpočet těchto vektorů jsem zvolil polynom stupně 12. To znamená, že pro výpočet jednoho vektoru bylo zapotřebí 13 časových okamžiků. Čas, pro který jsem vektor počítal se tedy nacházel uprostřed intervalu, který byl použit pro výpočet. Stavové vektory jsem počítal pro týdenní data. Tj. pro pozorování, která pokrývají dobu 7 dnů a to po 15 minutách. Vzhledem k tomu, že pro výpočet jednoho vektoru je třeba 13 časových okamžiků, nebylo možné určit vektory pro prvních a posledních šest časů pozorování. To znamená, že celkový počet stavových vektorů byl 7x24x4-6-6=660. Jako složku polohy stavového vektoru
jsem vzal souřadnice
družice v systému ICRS. Druhou složku vektoru, rychlost, jsem počítal jako derivaci polohy v daném souřadnici
zvlášť.
Konkrétní
způsob
čase a to pro každou jakým
jsem
vektory
20
vypočetl
je
následující.
Vektor
měření
(v
našem
případě
souřadnice družice):
x T = ( x1 , x2 , x3 ,..., x13 ) Vektor koeficientů funkce času je dán takto:
h T = (a0 , a1 , a2 ,..., a12 ) a konečně matice derivací funkce času je definována jako: ⎛1 ⎜ ⎜1 A = ⎜1 ⎜ ⎜M ⎜ ⎝1
Pro
větší
numerickou
t1
t12
t2
t 22
t3
t32
M
M
t13
t132
stabilitu
L t112 ⎞ ⎟ L t 212 ⎟ L t312 ⎟⎟ O M ⎟ 12 ⎟ L t13 ⎠
řešení
byl
výpočet
prováděn
v hodinách a kilometrech a zároveň sloupce 1-4 v matici A byly zvětšeny desetkrát. Vektor koeficientů h je tedy dán vztahem:
h = A −1 .x Pro souřadnice v daném časovém okamžiku (v tomto případě v časovém okamžiku t7) platí rovnice:
x = a0 + a1.t 7 + a2 .t 72 + ... + a12 .t 712 a pro rychlost družice v čase t7 je:
v=
[1.7]
dx = a1 + 2a2 .t 7 + 3a3 .t 72 + 4a4 .t 73 + ... + 12a12 .t 711 dt
Transformace souřadnic [x,y,z]ICRS družice a její rychlosti na Keplerovy dráhové elementy
Existuje odpovídají
šest šesti
Keplerových
dráhových
elementů,
integračním
konstantám,
tj.
které složkám
vektorového integrálu ploch a složkám vektorového integrálu Laplaceovu. Dělí se na vnější dráhové elementy – rektascenze
21 z
výstupního uzlu Ω, sklon roviny dráhy i, argument perigea ω -(obr. č.:9) a vnitřní
dráhové
poloosa
dráhy
elementy a,
pravá anomálie ν A
pro
dráhových
–
i
•
ρ
90o
N
v
hlavní
excentricita
ρ L
e,
90o
λ ω
- (obr. č.:10).
transformaci elementů
do
jsem
E
těchto
x
zvolil
i
Ω
y
I
obr. č.: 9
postup, při němž se používá výpočtu právě
η
zmíněných P’
B
P”
P
S a.e
v ω
b
a
E
P’’’
tento
postup
zdál,
Kromě co
do
toho
se
mi
programového
zajištění, nejvýhodnější.
M
Máme tedy šest počátečních podmínek.
ξ I
obr. č.: 10
vhodné
integrálů.
vypočítat
Jsou to souřadnice družice ICRS
(definují
⎛ dx j rychlosti ⎜⎜ ⎝ dt
pět
⎞ ⎟⎟ ⎠
integrálů.
(x ) j
její
polohu)
v epoše
t.
První
tři
v systému a
složky
Nejprve jsou
je
složky
vektorového integrálu ploch a jsou dány vektorovým součinem vektoru polohy a vektoru rychlosti:
()
⎛dx⎞ ⎟ L = x ⊗ ⎜⎜ ⎟ dt ⎠ ⎝
a v rozepsaném tvaru:
⎛ dx ⎞ ⎛ dx ⎞ L1 = ( x 2 ).⎜ 3 ⎟ − ( x3 ).⎜ 2 ⎟ ⎝ dt ⎠ ⎝ dt ⎠ ⎛ dx ⎞ ⎛ dx ⎞ L2 = (x3 ).⎜ 1 ⎟ − ( x1 ).⎜ 3 ⎟ ⎝ dt ⎠ ⎝ dt ⎠ ⎛ dx L3 = ( x1 ).⎜ 2 ⎝ dt
⎞ ⎛ dx ⎞ ⎟ − ( x 2 ).⎜ 1 ⎟ ⎠ ⎝ dt ⎠
Dále je třeba vypočítat integrál energie, hamiltonián H. Ten plyne, v případě konzervativních soustav, ze zákona zachování energie. H = T – V = konst.
22
kde T je kinetická energie družice o hmotnosti m:
T=
2 2 ⎤ 1 ⎡⎛ dρ ⎞ 2 ⎛ du ⎞ + m ⎢⎜ ρ ⎟ ⎜ ⎟ ⎥ 2 ⎢⎣⎝ dt ⎠ ⎝ dt ⎠ ⎥⎦
a V je její silová funkce:
V =m
GM
ρ
přičemž [ρ,u] jsou polární souřadnice družice v rovině její dráhy. Úhel u je úhel měřený od polární osy o (směřuje do výstupního uzlu I dráhy družice) ve směru pohybu družice. Vzhledem
k posledním
dvěma
uvedeným
vztahům
je
tedy
hamiltonián H roven: 2 2 ⎤ GM 1 ⎡⎛ dρ ⎞ 2 ⎛ du ⎞ H = m ⎢⎜ ⎟ + ρ ⎜ ⎟ ⎥ − m 2 ⎣⎢⎝ dt ⎠ ρ ⎝ dt ⎠ ⎦⎥
A
vzhledem
k tomu,
že
výraz
v hranaté
závorce
je
vlastně
výrazem pro rychlost družice, můžeme psát:
H=
GM 1 2 mv − m 2 ρ
a potom se bude hamiltonián H rovnat výrazu: 1 3 ⎛ dx j H = m∑ ⎜⎜ 2 j =1 ⎝ dt
Posledním
integrálem,
2
⎡ 3 ⎤ ⎞ ⎟⎟ − m.GM ⎢∑ (x j )2 ⎥ ⎠ ⎣ j =1 ⎦
který
je
třeba
−1 / 2
spočítat
je
vektorový
integrál Laplaceův. Ten je definován vektorovým vztahem: ⎡d ρ
⎤ ρ ⊗ L ⎥ − GM ρ ⎣ dt ⎦
λ=⎢
23
a rozepsán ve složkách má tvar:
dx dx 2 − L2 3 − GM dt dt dx dx λ 2 = L1 3 − L3 1 − GM dt dt dx dx λ3 = L2 1 − L1 2 − GM dt dt
x1
λ1 = L3
ρ
x2
ρ x3
ρ
A nyní můžeme přejít rovnou na výpočet Keplerových dráhových elementů. Z vektorového integrálu ploch dostáváme hned několik výrazů pro výpočet rektascenze výstupního uzlu Ω a sklonu i:
(L tan i =
2 1
L tan Ω = − 1 L2 sin Ω. sin i =
L1 L
− cos Ω. sin i =
+ L22 L3
)
L2 L
cos i =
L3 L
Excentricitu můžeme spočítat podle vzorce:
e=
λ GM
Velká poloosa je dána vztahem:
L2 1 GM a=− m nebo také a = 2 H GM 1 − e 2
(
)
Pravou anomálii můžeme spočítat dle vzorce: 1⎡ 3 2⎤ cos v = ⎢∑ (x j ) ⎥ λ ⎣ j =1 ⎦
−1 / 2
3
∑ λ .x j =1
j
j
Posledním Keplerovým dráhovým elementem je argument perigea. Ten určíme na základě skalárního součinu Laplaceova vektoru λ a vektoru uzlové přímky (cos Ω; sin Ω; 0): cos ω =
λ1 λ cos Ω + 2 sin Ω λ λ
24
Vzhledem k tomu, že dráhy družic GPS jsou téměř kruhové, bude výhodnější samostatně
interpolovat argument
argument
perigea
a
u=ω+ν
deklinace
pravou
anomálii.
než
Argument
deklinace je dán vztahem:
⎡1 ⎤ cos u = ⎢ ( x1 . cos Ω + x 2 . sin Ω )⎥ ⎣r ⎦ Další
úprava
interpolovat
pro hlavní
výpočty poloosu
spočívá dráhy
v tom,
družic
(a
že tedy
nebudeme i
jejich
excentricitu), ale přímo průvodič družic: r = a.(1 − e. cos E )
přičemž excentrickou anomálii E spočteme ze vztahu:
tan E =
[1.8]
(x • v ).
2.GM − v2 x x⎞ ⎛ GM .⎜1 − ⎟ ⎝ a⎠
Interpolace terestrických souřadnic [x,y,z]ITRS
Vzhledem k časovému průběhu souřadnic (viz. obr. č.: 11a, 11b a 11c) jsem se rozhodl nepoužít polynomickou funkci, ale interpolaci pomocí Fourierovy harmonické řady. Ta je popsána
25000
P rů b ě h so u řad n ice X v čase (G P S týd e n 899, d ru ž ice č.: 1)
15000 10000 5000 0 -1 0 0 0 0 -1 5 0 0 0 -2 0 0 0 0 -2 5 0 0 0 č a s (p o č á te k 5 0 5 3 7 .0 0 0 0 MJ D = 1 , 1 d íle k = 1 5 m in ) o b r. č .: 1 1 a
600
500
400
300
200
100
-5 0 0 0
0
km (vztaženo ke geocentru)
20000
25
v kapitole [1.2]. Z grafů (obr. č.: 11a, 11b) je patrné, že časová závislost dat,
v tomto
periodicky
případě
se
souřadnice
opakujících
X
funkcí.
a
Y,
Dále
je
součtem
můžeme
dvou
odhadnout
periody obou funkcí (Ph je hlavní perioda a Pv je vedlejší Průběh souřadnice Y v čase (GPS týden 899, družice č.: 1) 25000 15000 10000 5000 0 600
500
400
300
200
100
-5000
0
km (vztaženo ke geocentru)
20000
-10000 -15000 -20000 -25000 čas (počátek 50537.0000 MJD = 1, 1 dílek = 15 m in) obr. č.: 11b
perioda) a hlavní amplitudu. Vedlejší amplituda se z těchto grafů odhaduje velmi spatně. Číselné hodnoty těchto veličin jsou
uvedeny
v kapitole
pojednávající
o
praktickém
řešení
jednotlivých problémů. Tyto hodnoty nejsou pouze odhadem, ale jsou podpořeny i analýzou pomocí Fourierovy transformace. Avšak
z informací,
interpolační
funkci
které
(pro
nyní
máme,
interpolaci
můžeme
sestavit
v souřadnicích
X,Y),
která by se co nejlépe přimykala danému časovému průběhu dat. Je to funkce:
(
)
(
y = a.sin t (1) + A + b.sin t (2 ) + B
)
kde neznámé jsou: a jako hlavní perioda, A je hlavní fázový posun, b je vedlejší perioda a konečně B je vedlejší fázový posun. Dále, jak je popsáno v kapitole [1.2], jsou t(1) a t(2) (čísla v závorkách nejsou mocniny, ale indexy odlišující od sebe obě proměnné) nové proměnné, které dostaneme z původních proměnných x takto:
26
t i(1) =
2π 2π xi a ti(2 ) = xi Pv Ph
Využijeme nyní součtový vzorec pro součet argumentů ve funkci sinus a dále zavedeme substituci:
a.sin A = A1 a. cos A = A2
,
b. sin B = B1 b. cos B = B2
Potom bude konečný tvar interpolační funkce tento:
y = A1 . cos t (1) + A2 .sin t (1) + B1 . cos t (2 ) + B2 .sin t (2 ) Je zřejmé, že funkce pro interpolaci v souřadnicích Z (obr. č.:
11c)
bude
jednodušší.
Bude
obsahovat
pouze
první
člen
Průbě h souřadnice Z v čase (G PS týde n 899, družice č.: 1) 25000 15000 10000 5000 0 600
500
400
300
200
100
-5000
0
km (vztaženo ke geocentru)
20000
-10000 -15000 -20000 -25000 čas (počáte k 50537.0000 MJD = 1, 1 dílek = 15 m in) obr. č.: 11 c
Fourierovy
harmonické
řady.
Její
odvození
zde
ale
uvádět
nebudu, neboť je naprosto analogické s výše odvozenou funkcí pro souřadnice X a Y, s omezením na jeden člen řady. Konečný tvar tedy bude: y = A1 . cos t + A2 . sin t
kde
ti =
2π xi P
27
[1.9]
Interpolace Keplerových dráhových elementů
Interpolační funkce dráhových elementů jsou součtem dvou dílčích
funkcí.
případech
i
Jsou
to
funkce
konstantní),
polynomické
které
vystihují
elementů po dobu, kterou zabírají týden),
a
pak
jsou
to
funkce
(v
některých
trend
dráhových
analyzovaná data (jeden GPS harmonické,
které
vystihují
harmonický časový průběh dráhových elementů samotných. V některých případech je použita funkce polynomická pro vystižení
týdenního
trendu,
ale
ve
skutečnosti,
vezmeme-li
v potaz data většího časového úseku než jeden týden, přejde tento trend na funkci harmonickou, jejíž perioda je mnohem větší
než
právě
jeden
týden.
A
dlouhodobý
časový
vývoj
dráhových elementů by pak měl ještě další průběhy. Konečné
podoby
interpolačních
funkcí
uvedené
v této
kapitole jsou výsledkem několikadenních pokusů o co nejlepší vystižení použití [1.9]
průběhu analýzy
a
pomocí
v jejích
interpolačních grafů
časové
Fourierovy
podkapitolách
funkcí,
časových
závislosti
které
průběhů
dráhových
transformace.
budu
budou
elementů
uvádět
jednotlivých
V kapitole
konečné
komentovány dráhových
za
na
tvary
základě elementů.
Komentáře založené na analýze pomocí Fourierovy transformace budou uvedeny dále v textu v kapitole [2.3], která se týká praktických
výpočtů
a
mezivýsledků.
Tam
také
budou
uvedeny
číselné hodnoty period a amplitud jednotlivých harmonických průběhů. Ještě je třeba podotknout, že pro praktickou část této práce,
tedy
pro
výpočty,
jsem
si
náhodným
výběrem
zvolil
družici PRN1. Nicméně vše co je v této práci uvedeno platí i pro všechny ostatní družice GPS. Výpočty byly prováděny ještě pro některé další družice. Výsledky výpočtů ale nemohu uvádět a to z toho důvodu, že by tato práce nabyla ohromných rozměrů.
28
Rektascenze výstupního uzlu dráhy družice - Ω
[1.9.1]
Z grafu časového průběhu Ω (obr. č.: 12) je zřejmé, že interpolační funkce bude mít tvar součtu funkce lineární a harmonické. potvrdila,
Dále že
je
zde
vidět,
figuruje
a
Fourierova
pouze
[ o]
Rektascenze výstupního uzlu dráhy Ω
jedna
transformace harmonická
funkce
DETAIL
158,250
158,325 158,300 158,275 158,250 158,225 158,200 158,175 158,150 158,125 158,100 158,075 158,050
to
158,225
600
500
400
300
200
100
0
158,200
(s periodou P) a že tedy ve výsledné interpolační funkci bude jen jeden harmonický člen. Interpolační funkce tedy bude vypadat takto: y = a.x + b + c.sin (t + C )
Pomocí součtového vzorce pro součet argumentů funkce sinus a substituce:
c. sin C = C1 c. cos C = C2 a
po
příslušných
úpravách
dostáváme
interpolační funkce: y = a.x + b + C1 . cos t + C 2 .sin t
přičemž t bude rovno výrazu:
ti =
2π xi P
konečnou
podobu
320
300
280
260
240
220
200
čas [počátek 50537.0000 MJD=1, 1 dílek po 15 minutách] obr. č.: 12
29
[o]
[1.9.2]
Sklon dráhy družice - i
S k lo n d rá h y i
54,706 54,705 54,704 54,703 54,702 54,701 54,700 54,699 54,698 54,697 54,696
600
500
400
300
200
100
0
č a s [p o č á t e k 5 0 5 3 7 . 0 0 0 0 M JD = 1 , 1 d íle k p o 1 5 m in u tá c h ] o b r. č .: 1 3
Jak je vidět z týdenních dat (obr. č.:13) má sklon dráhy družice
průběh
přiblížila
harmonický.
tomuto
harmonických
průběhu
funkcí
o
Funkce, by
která
vypadala
různých
by
jako
periodách.
se
nejlépe
součet
Nicméně
dvou hlavní
perioda tohoto časového průběhu je delší než jeden týden a proto funkce, kterou jsem nakonec zvolil jako interpolační je součtem
kubické
periodou
P)
polynomické
s jedním
funkce
harmonickým
a
harmonické členem.
funkce
Konečný
(s
tvar
interpolační funkce tedy je:
y = a.x 3 + b.x 2 + c.x + d + e.sin (t + E ) Součtovým
vzorcem
pro
součet
argumentů
zavedením substituce:
e.sin E = E1 e. cos E = E 2 získáme finální podobu funkce:
y = a.x 3 + b.x 2 + c.x + d + E1 . cos t + E 2 .sin t přičemž t bude opět rovno výrazu:
ti =
2π xi P
funkce
sinus
a
30
[1.9.3]
Průvodič družice - r
P rů v o d ič d rá h y r
26680 26650 26620 [km]
26590 26560 26530 26500 26470 26440 600
500
400
300
200
100
0
č as [poč átek 50537.0000 M JD = 1, 1 dílek po 15 m inutác h] obr. č .: 14
Již
z grafu
časového
průběhu
velikosti
průvodiče
dráhy
družice (obr. č.:14) je vidět, že interpolační funkce bude velmi jednoduchá. Trend tohoto časového průběhu průvodiče je konstantní. Je zřejmé, že v harmonickém průběhu figuruje pouze jedna frekvence s periodou Ph a amplitudou b. Výsledná podoba interpolační funkce tedy je: y = a + b. sin (t + B )
kde
b je zmiňovaná hlavní amplituda,
B je patřičný fázový
posun a t je nová proměnná získaná transformací:
ti =
2π xi Ph
Za využití již výše několikrát zmíněného součtového vzorce a substituce:
b. sin B = B1 b. cos B = B2 dostaneme
po
jednoduchých
úpravách
interpolační
tvaru: y = a + B1 . cos t + B2 . sin t
funkci
ve
31
[1.9.4]
Argument deklinace dráhy družice - u A rg u m e n t d e k lin a c e u
180 150 120 o
[]
90 60 30 0 600
500
400
300
200
100
0
čas [počátek 50537.0000 M JD = 1, 1 dílek po 15 m inutách] obr. č.: 15
Z grafu (obr. č.:15) je patrné, že veličinou změny
(stejně
své
velikosti
nejostřejší. zmíněném
jako
průvodič),
v čase.
Závislosti,
grafu,
nejlépe
A
kterou
argument deklinace je
která nejen
můžeme
tvarově
vykazuje
největší
největší
ale
také
ve
výše
pozorovat
vyhovuje
pilovitá
funkce,
která je dána součtem řady:
y =a+
8.b ⎧ 1 1 ⎫ sin[t + B ] − 2 sin[3(t + B )] + 2 sin[5(t + B )] − ...⎬ 2 ⎨ 3 5 π ⎩ ⎭
z níž pro náš účel bohatě postačí první tři členy této řady. Koeficient
a
vyjadřuje
konstantní
trend
časové
závislosti,
koeficient b je pak její amplituda a B má význam fázového posunu.
Proměnná
t
je
pak
opět
nová
proměnná
vzniklá
transformací, která má tvar:
ti = Následují
analogické
v předcházejících
2π xi Ph
substituce
podkapitolách
a
analogické
kapitoly
úpravy
[1.9],
po
získáme konečnou podobu interpolační funkce: y = a + B2 . sin t + B1 . cos t − C 2 . sin 3t − C1 cos 3t + D2 sin 5t + D1 cos 5t
jako nichž
32
[2] PRAKTICKÁ REALIZACE VÝPOČTŮ V této konečné
kapitole
výsledky
a
popíši
postupy
okomentuji
výpočtů,
časové
uvedu
závislosti
na
dílčí
i
základě
analýzy dat pomocí Fourierovy transformace. Závěr a komentáře výsledků pak budou uvedeny v kapitole [3]. Dále
zde
budou
uvedeny
ukázky
datových
souborů
se
souřadnicemi družic a souřadnicemi pólu a jinými důležitými parametry
vstupujícími
do
výpočtů.
Popíši
zde
také
způsob,
kterým jsem získal data potřebná pro výpočetní práce. Pro úplnost bych chtěl podotknout, že, kromě Fourierových transformací (ty jsem počítal v programu Matlab 4.0), jsem vše počítal na software, který jsem si sám naprogramoval v jazyce C++ za použití vývojového prostředí C++ Builder v3.0. Zdrojové kódy, programové soubory (*.exe) a další dokumenty týkající se této práce jsou přiloženy na CD. Popis programů je v příloze č.:3.
[2.1]
Vstupní datové soubory a jejich struktura
Datové
soubory,
které
obsahují
potřebné
informace
jsem
získal pomocí Internetu a služby FTP (File Transfer Protocol) a to z vládního serveru organizace NASA (National Aeronautic and
Space
Administration
kosmický
prostor)
cddis.gsfc.nasa.gov
–
Národní
v USA. a
adresář,
správa
Adresa kde
se
pro
letectví
serveru data
a
je:
nacházejí
je:
gps3:[products]/“week“, kde „week“ je číslo týdne (například: 0899), který nás zajímá. V těchto
adresářích
se
tedy
dají
najít
jak
soubory
s predikovanými, tak i s finálními souřadnicemi družic a také predikované a finální souřadnice pravého pólu v soustavě CIO (viz. kapitola [1.5.4]). Soubory IGPxxxxy.sp3 jmenují
s predikovanými a
soubory
IGSxxxxy.sp3
souřadnicemi
s finálními
družic
souřadnicemi
se
jmenují
družic
se
a dále soubor s finálními souřadnicemi
pravého pólu má jméno IGSxxxx7.erp a soubory s predikovanými
33
souřadnicemi
pravého
pólu
se
jmenují
IGPxxxxy.erp.
Písmena
v názvech souborů mají tento význam: xxxx
-
číslo GPS týdne (například 0899)
y
-
číslo dne v týdnu od 0=neděle až do 6=sobota
Každý ze souborů má svou „hlavičku“ a strukturované „tělo“. Příklad „hlavičky“ souboru souřadnic družic je následující: #aP1997 3 30 0 0 .00000000 96 ORBIT ITR94 HLM IGS ## 899 .00000000 900.00000000 50537 .0000000000000 + 25 1 2 3 4 5 6 7 9 10 14 15 16 17 18 19 21 22 + 23 24 25 26 27 29 30 31 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ++ 4 5 5 5 4 5 5 5 4 5 5 5 5 5 5 5 4 ++ 5 4 5 4 5 5 4 5 0 0 0 0 0 0 0 0 0 ++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 %c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc %c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc %f .0000000 .000000000 .00000000000 .000000000000000 %f .0000000 .000000000 .00000000000 .000000000000000 %i 0 0 0 0 0 0 0 0 0 %i 0 0 0 0 0 0 0 0 0 /* FINAL ORBIT COMBINATION FROM WEIGHTED AVERAGE OF: /* cod emr esa gfz jpl ngs sio /* REFERENCED TO GPS CLOCK AND TO WEIGHTED MEAN POLE: /*
Údaje které byly pro mne důležité a které lze vyčíst z této „hlavičky“ jsou: datum:
30.3.1997
souřadný systém:
ITRF94, zpracováno HLM IGS
čas:
.00000 sekunda týdne 899 (neděle)
Dalšími údaji jsou počet zpracovávaných družic (25), jejich PRN (Pseudo Random Nubers - pseudonáhodná čísla), MJD týkající se půlnoci daného dne, tedy 30.3.1997 0h 0m 0s. Tělo souboru se souřadnicemi družic má pevně stanovenou strukturu. Ta je velmi dobře patrná * P P P P P P
1997 3 30 0 0 .0000 1 -20672.111626 15834.955147 5355.398110 2 11860.349918 -15541.883699 -17345.088179 3 -4894.172170 18383.557550 -18636.090602 4 8962.781399 -12072.258233 22011.471799 5 -12508.920324 -16608.266153 16568.397184 6 -22807.783717 1606.657070 13838.093491
z ukázky:
16.308960 -349.195699 113.435487 20.183024 89.704742 4.361641
34
Jsou
zde
uvedeny
PRN
0h
0s
(30.3.1997,
0m
v prvním
řádku
obsahuje
souřadnice
korekci
jejich
jednotlivých =
každého
50537.0000 bloku
družic,
hodin
v
družic
a
je
které
pro
MJD),
který
označen jsou
daný
mikrosekundách.
je
„*“,
seřazeny Bloky
okamžik uveden
dále dle
jsou
blok PRN
a
seřazeny
v časové souslednosti vždy po 15 minutách. Datový
soubor
s finálními
souřadnicemi
pravého
pólu
a
korekcí (UT1-UTC) má tuto strukturu: Source: Xpole,Ypole,Xrt,Yrt: weighted average of centres; LODsig == 0 -> LOD,UT1-UTC: from IERS Bulletin A ; LODsig != 0 -> LOD: weighted average of centres; UT1-UTC: integrated from the most current (non predicted) Bull. A value. Orbits: to be used with the IGS Final Orbits MJD Xpole Ypole UT1-UTC LOD Xsig Ysig UTsig LODsig Nr Nf Nt Xrt Yrt Xrtsig Yrtsig (10**-5") usec usec (10**-5") usec usec (10**-5"/d) (10**-5"/d) 50537.50 -19091 32491 -287709 2395 4 6 47 14 0 0 0 -12 314 4 13 50538.50 -19102 32796 -290161 2397 8 11 42 18 0 0 0 -35 265 26 32 50539.50 -19067 33139 -292528 2518 5 13 43 19 0 0 0 76 319 13 15 50540.50 -18980 33523 -295143 2612 8 13 39 8 0 0 0 70 416 26 13 50541.50 -18875 33907 -297885 2710 8 14 45 13 0 0 0 102 376 11 12 50542.50 -18750 34290 -300588 2883 5 11 44 10 0 0 0 121 337 10 14 50543.50 -18655 34635 -303502 2974 9 23 46 15 0 0 0 105 353 9 15
Informace,
které
byly
pro
řešení
daných
problémů
potřebné
jsou: MJD a pro něj uváděné souřadnice pólu v soustavě CIO a korekce
hodin
(UT1-UTC).
Z dalších
informací
obsažených
v tomto souboru uvedu například rozdíl aktuální délky dne a 24 hodin v sekundách, charakteristiky přesnosti souřadnic pólu, rozdílu (UT1-UTC), délky dne (LOD).
[2.2]
Interpolace terestrických souřadnic [x,y,z]ITRS
V kapitole transformaci frekvencí,
[1.8]
jako
period
na a
jsem
se
nástroj
amplitud
odvolával mnohem
na
Fourierovu
přesnějšího
harmonických
časových
určení průběhů
dat. Fourierovy transformace tedy pro jednotlivé terestrické souřadnice vypadají takto:
35
DETAIL 1
1
0
0
-1
amplituda [*106]
amplituda [*106]
Fourierova transformace souřadnice X (GPS týden 899, družice č.:1)
-1 -2 -3 -4 -5
-2 -3 -4 -5
-6
30
5 4 amplituda [*106]
4 amplituda [*106]
30
25
DETAIL
5
3 2 1 0
3 2 1 0
-1
frekvence
DETAIL
amplituda [*105]
9 8 7 6 5 4 3 2 1 0 -1
9 8 7 6 5 4 3 2 1 0 -1 35
30
frekvence
25
20
15
10
5
0
600
500
400
300
200
100
0
frekvence
15
Fourierova transformace souřadnice Z (GPS týden 899, družice č.: 1)
10
5
0
600
500
400
300
200
100
-1
frekvence
amplituda [*105]
25
20
frekvence
Fourierova transformace souřadnice Y (GPS týden 899, družice č.: 1)
0
20
15
frekvence
10
5
0
600
500
400
300
200
100
0
-6
36
Informace, které z těchto grafů (a z grafů na obr. č.:11a, 11b, 11c) můžeme vyčíst je možné sestavit do tabulky: souřadnice
X
Y
Z
hlavní
7
7
14
vedlejší
21
21
---
hlavní [měř.]
96
96
48
vedlejší [měř.]
32
32
---
20880
21600
5540
---
přibližné
fázové
frekvence perioda amplituda
hlavní [km] 20880 vedlejší [km]
5540
tab. č.: 1
Hodnoty
v tabulce
a
odečtené
z grafů na obr. č.:11a, 11b, 11c přibližných
koeficientů
posuny
byly použity pro výpočet
aproximačních
funkcí,
které
jsou
uvedeny v kapitole [1.8]. Z grafů Fourierových transformací jednotlivých souřadnic družic je dále vidět to, že dochází k velmi ostrým změnám hodnot souřadnic. To je patrné z toho, že píky ostatní signál (na časový průběh souřadnic se můžeme dívat jako na určitý signál)
výrazně
převyšují.
Přechody
mezi
píky
a
ostatním
signálem jsou velmi ostré a úzké. To znamená, že časový průběh souřadnic není zatížen náhodnými rušivými vlivy, tzv. šumem, a je tedy z pohledu signálu poměrně čistý. Matice
plánu
interpolačních
byla
funkcí
sestavena dle
z parciálních
jednotlivých
derivací
neznámých,
tedy
koeficientů interpolačních funkcí. Vektor neznámých pak tvoří jednotlivé koeficienty těchto funkcí (viz. aplikace kapitoly [1.2]).
Výsledné
funkcí
a
odhady
vektory jejich
v tabulkách č.:2a, 2b, 2c.
neznámých
parametrů
středních
chyb
interpolačních jsou
sestaveny
37
souřadnice X neznámé
hi [km]
mhi [km]
A1
-15280,46
255,90
A2
-14235,64
255,90
B1
-5360,39
255,90
B2
1414,43
255,90
tab. č: 2a
souřadnice Y neznámé
hi [km]
mhi [km]
A1
14437,66
258,49
A2
-15274,42
258,49
B1
1420,67
258,49
B2
5451,62
258,49
tab. č: 2b
souřadnice Y neznámé
hi [km]
mhi [km]
A1
5335,97
422,51
A2
-20925,84
422,51
tab. č: 2c
Z hodnot v tabulkách 2a, 2b a 2c jsem dále vypočítal původní koeficienty interpolačních funkcí pomocí umocnění a sečtení resp.
podělení
substitučních
výrazů,
které
jsou
popsány
v kapitole [1.8]. Číselné hodnoty koeficientů jsou přehledně uspořádány do tabulek: souřadnice X neznámé
hi
mhi
a [km]
20884,106
255,900
A [rad]
3,962374351
0,012253337
b [km]
5543,861
255,900
B[rad]
4,970375753
0,046159167
38
souřadnice Y neznámé
hi
mhi
a [km]
21017,943
258,490
A[rad]
2,384349331
0,012298539
b [km]
5633,690
258,490
B[rad]
0,2549261762
0,045882893
souřadnice Y
[2.3]
neznámé
hi
mhi
a [km]
21595,448
422,510
A [rad]
2,891919021
0,019564772
Interpolace Keplerových dráhových elementů
K analýze
Keplerových
dráhových
elementů
opět
použiji
Fourierovu transformaci. Na základě interpretace grafů této transformace a elementů
grafů časových průběhů jednotlivých dráhových
určíme
interpolačních
přibližné
funkcí
a
hodnoty
kromě
jednotlivých
toho
nám
též
koeficientů
pomůže
zjistit
teoretickou použitelnost této metody pro určení predikovaných souřadnic
družic
na
základě
určení
interpolačních
funkcí
Keplerových dráhových elementů a jejich koeficientů. Kapitola je rozdělena do jednotlivých podkapitol, v nichž jsou popsány grafy Fourierových analýz Keplerových dráhových elementů. Dále jsou v podkapitolách uvedeny frekvence, periody a
amplitudy
Keplerových
dráhových
elementů.
V poslední
podkapitole jsou pak uvedeny tabulky s výslednými hodnotami koeficientů
jednotlivých
středními chybami.
interpolačních
funkcí
a
jejich
39
Rektascenze výstupního uzlu dráhy družice - Ω
[2.3.1]
Graf Fourierovy transformace časového průběhu rektascenze výstupního uzlu dráhy družice má tuto podobu: Fourierova transformace - Ω (GPS týden 899, družice č.: 1)
0,40
0,35
0,35
0,30
0,30
0,25
0,25
amplituda
amplituda
0,40
DETAIL
0,20
0,20
0,15
0,15
0,10
0,10
0,05
0,05
0,00
Z tohoto grafu a z grafu na obr. č.:12 je možné vyčíst tyto důležité informace:
Ω
Keplerův dráhový element frekvence
hlavní
1
vedlejší
27
hlavní [měř.]
---
vedlejší [měř.]
24
hlavní [o]
---
perioda amplituda
vedlejší [o]
0,0016
tab. č.: 3
Vzhledem (opět
se
na
k tomu,
že
píky
časový
průběh
nad
ostatním
Keplerových
pozadím
dráhových
signálu elementů
můžeme dívat jako na jistý signál a danou periodou, frekvencí a amplitudou) nejsou příliš převýšeny, můžeme konstatovat, že časová změna hodnot rektascenze výstupního uzlu dráhy družice nenabývá závratných velikostí. Dále je možné říci, že časový
50
frekvence
45
40
35
30
25
20
15
frekvence
10
5
0
600
500
400
300
200
100
0
0,00
40
průběh této veličiny není zatížen šumem a je tedy poměrně čistý. To je patrné z toho, že přechody do píků jsou ostré. Rád bych ještě okomentoval tu skutečnost, že Fourierova transformace odhalila ještě jednu frekvenci, jejíž perioda je vzhledem k týdennímu časovému intervalu dat dlouhodobá. Jedná se s největší pravděpodobností a kombinaci jevů, které mají periodu cca. 14 a 28 dní a amplitudu v součtu okolo 9“. Tyto jevy
se
ale
v mých
výpočtech
neuplatní,
neboť
mají
na
rektascenzi výstupního uzlu v horizontu jednoho týdne nepatrný vliv.
[2.3.2]
Sklon dráhy družice - i
Na obrázku č:13 můžeme vidět časový průběh sklonu dráhy družice, jehož graf Fourierovy transformace vypadá takto: DETAIL 0,4
0,4 0,3 0,2 0,1 0 -0,1 -0,2 -0,3 -0,4 -0,5
0,3 0,2 amplituda
amplituda
Fourierova transformace - i (GPS týden 899, družice č.: 1)
0,1 0 -0,1 -0,2 -0,3 -0,4
uvedu
tabulku
(tab.
č.:4)
s hodnotami
frekvencí, period a amplitud jevů, které odhalila Fourierova analýza a graf časového průběhu sklonu dráhy družice:
50
frekvence
45
40
35
30
25
20
nejprve
15
Opět
10
frekvence
5
0
600
500
400
300
200
100
0
-0,5
41
Keplerův dráhový element frekvence
i
hlavní
1
vedlejší
27
hlavní [měř.]
---
vedlejší [měř.]
24
hlavní [o]
---
perioda amplituda
vedlejší [o]
0,0013
tab. č.: 4
Fourierova
transformace
odhalila
i
zde
jev,
který
má
frekvenci 1 (má delší periodu než týden). Tento jev jsem již popsal v kapitole [1.9.2] a je to právě ten trend, který jsem nahradil
kubickým
polynomem.
Vzhledem
k týdennímu
časovému
horizontu dat je tento trend dlouhodobější. Nicméně s ním bylo nutno počítat, protože jeho amplituda má velikost 9” (perioda činí 14 dní). Ani
v tomto
případě
píky
podstatně
nepřevyšují
ostatní
pozadí signálu. To znamená, že amplituda harmonického průběhu hodnot
této
veličiny
nebude
velká
a
že
nebude
docházet
k dramatickým změnám hodnot této veličiny v čase. Nicméně je vidět, že průběh sklonu dráhy družice v čase je mírně zatížen náhodnými
rušivými
vlivy,
tzv.
šumem,
protože
přechod
mezi
píkem a pozadím signálu je poměrně pozvolný. [2.3.3]
Průvodič družice - r
Fourierova grafu
na
následující
transformace hodnoty
transformace
a
za
stránce.
pomoci
frekvence,
této
veličiny
Na
obrázku
periody
a
základě č.:14
je
této
jsem
amplitudy
znázorněna
Keplerův dráhový element frekvence perioda amplituda
časového r
hlavní
14
hlavní [měř.]
48
hlavní [km]
Fourierovy
stanovil
průvodiče:
99,408
na
tyto
průběhu
42
DETAIL
Fourierova transformace - r (GPS týden 899, družice č.: 1) 45 amplituda [*105]
amplituda [*105]
35 25 15 5 -5 -15
50
45
40
35
transformace
48
měření
z
grafu
Dále
30
99,408 km.
frekvence
Fourierovy
s periodou
25
průběh
20
Z grafu
15
frekvence
10
5
0
600
500
400
300
200
100
0
40 35 30 25 20 15 10 5 0 -5 -10 -15
(po
lze
15
vyčíst
minutách)
Fourierovy
periodický
a
amplitudou
transformace
můžeme
usoudit na to, že dochází k velkým změnám hodnot této veličiny v čase,
protože
pík
velmi
převyšuje
pozadí
signálu.
Dále
přechod mezi píkem a pozadím signálu není ostrý, ale poměrně pozvolný.
Můžeme
tedy
říci,
že
časový
průběh
hodnot
této
veličiny je mírně zatížen šumem. [2.3.4]
Argument deklinace dráhy družice - u
Graf Fourierovy transformace je vyobrazen na obrázku: DETAIL
Fourierova transformace - u (GPS týden 899, družice č.: 1)
2,5
2,5
2,0 amplituda [*104]
amplituda [*104]
2,0 1,5 1,0 0,5 0,0
1,5 1,0 0,5 0,0
-0,5
-0,5
-1,0
-1,0 100
90
80
frekvence
70
60
50
40
30
20
10
0
600
500
400
300
200
100
0
frekvence
43
Z grafu Fourierovy transformace argumentu deklinace je na první pohled zřejmé, že zde dochází k velmi rychlým a velkým změnám hodnot této veličiny v čase, protože zde píky dosahují značných hodnot. Dále vidíme, že přechod mezi píkem a zbylým pozadím
není
ostrý
a
proto
můžeme
konstatovat,
že
hladký
průběh bude mírně zatížen šumem. Tento šum ovšem nedosahuje takových
hodnot,
aby
velkou
měrou
ovlivňoval
výsledky
interpolace této veličiny. Je zde evidentní jedna frekvence o hodnotě
14,
jejíž
amplituda
mnohonásobně
převyšuje
ostatní
píky, které jsou mimo jiné také důkazem přítomnosti šumu, o kterém jsem se zmínil výše v tomto odstavci a jímž je zatíženo pozadí signálu. Hodnoty, které jsem vyčetl jak z Fourierovy transformace tak i z grafu časového průběhu pravé anomálie jsou uvedeny v následující tabulce: Keplerův dráhový element frekvence perioda
hlavní
14
hlavní [měř.]
48
o
amplituda
[2.3.5]
hlavní [ ]
89,7
Výsledky výpočtů
V této vepsány
u
kapitole
výsledky
budou
výpočtů
uvedeny
a
koeficientů
přehledně
do
tabulek
interpolačních
funkcí,
které jsou popsány v kapitole [1.9]. Pro stanovení přibližných hodnot těchto koeficientů jsem použil nejen hodnoty uvedené v kapitole [2.3], ale také odhady fázových posunů jednotlivých časových
závislostí.
Tyto
posuny
jsou
velmi
dobře
patrné
z grafů časových závislostí (obrázky č.:12-15) a proto jsem je do tabulek neuváděl. Tam kde tyto posuny patrné přímo nebyly, byl jsem odkázán pouze na jejich odhad. Nicméně do výpočtu přibližných koeficientů jsou tyto posuny zahrnuty.
44
Výsledky, ke kterým jsem při výpočtech došel tedy jsou: Keplerovy dráhové elementy
Ω
i
h [rad]
h [rad]
mh [rad]
a= -7,234447327E-6 1,865367182E-9 2,763452744 7,116319521E-7 b= C1= 2,371446143E-5 5,026123326E-7 C2= 1,384354289E-5 5,026777371E-7
mh [rad]
a= -1,61863788E-12 1,320797908E-9 b= c= -8,295882344E-8 0,9546683787 d= E1= 1,098019246E-5 E2= -1,929642108E-5
1,874316E-14 1,884701061E-11 5,363703836E-9 4,097804064E-7 1,440156275E-7 1,441275214E-7
Keplerovy dráhové elementy r
u
h [m] a= 26561121,466349 -50754,690849 B1= -84941,615005 B2=
Dále
h [rad]
mh [m]
jsem
191,188433 270,572742 270,153519
spočetl
mh [rad]
1,5731639219 a= B2= -0,6484541852 B1= 1,0919048630 C2= 0,1375101445 C1= 0,0047597538 D2= -0,0205526854 D1= -0,0435551059
původní
koeficienty
0,0028268662 0,0039945568 0,0040007570 0,0039941617 0,0040001698 0,0039938277 0,0039999081
v interpolačních
funkcích z nichž jsem vycházel v kapitolách [1.9.1] až [1.9.4] a to umocněním a sečtením, popřípadě podělením, jednotlivých substitučních jejich
výrazů.
středních
Výsledné
chyb
jsou
hodnoty
koeficientů
sestaveny
s odhady
v následujících
tabulkách: Keplerovy dráhové elementy
Ω h [rad] a= b= c= C=
i
-7,234447327E-6 1,865367182E-9 2,763452744
7,116319521E-7
0,000027459413 0,00000050262 1,042403945
h [rad]
mh [rad]
0,01830560308
a= b= c= d= e= E=
-1,61863788E-12
mh [rad] 1,874316E-14
1,320797908E-9
1,884701061E-11
-8,295882344E-8
5,363703836E-9
0,9546683787
4,097804064E-7
0,0000222017228
0,0000001441002
2,624258519818
0,006487920134
45
Keplerovy dráhové elementy r
u
h
mh
26561121,466349
a[m]= 98950,071766 b[m]= B[rad]= 4,1737917338
191,188433 270,263878 0,0016892822
a= b= B=
h [rad]
mh [rad]
1,5731639219
0,0028268662
1,5667263831
0,0032415819
2,1066994269
0,0046009995
[3] ZÁVĚR: Veškeré číselné hodnoty výpočtů koeficientů interpolačních funkcí jsou uvedeny v příslušných tabulkách v kapitolách a
[2.3.5].
V této
kapitole
bych
rád
okomentoval
[2.2]
praktickou
použitelnost těchto dvou metod, kterými jsem se v diplomové práci zabýval. Výsledky výpočtů uvedené v kapitole [2.2] ukazují přímo střední chyby s jakými by byly určeny souřadnice družice PRN1 z vypočítaných vyrovnaných koeficientů interpolačních funkcí. Jejich velikosti jsou: mx = 361,9 km my = 365,6 km mz = 422,5 km Pro zjištění odhadů středních chyb souřadnic družice PRN1 interpolovaných na základě Keplerových dráhových elementů je třeba
malého
zamyšlení.
Transformace
Keplerových
dráhových
elementů na kartézské souřadnice v systému ICRS je dána těmito vzorci:
x = r.(cos u. cos Ω − sin u. sin Ω. cos i ) y = r.(cos u. sin Ω + sin u. cos Ω. cos i ) z = r. sin u. sin i Podívejme se jak budou vypadat parciální derivace souřadnice x (u ostatních souřadnic jde o analogii a proto se budu v těchto
46
rozborech
zabývat
pouze
souřadnicí
x)
podle
jednotlivých
veličin: ∂x = (cos u. cos Ω − sin u. sin Ω. cos i ) ∂r ∂x = r.(− sin u. cos Ω − cos u. sin Ω. cos i ) ∂u ∂x = r.(− cos u. sin Ω − sin u. cos Ω. cos i ) ∂Ω ∂x = r.(sin u. sin Ω. sin i ) ∂i
Podle zákona hromadění středních chyb platí: 2
2
2
2
⎛ ∂x ⎞ ⎛ ∂x ⎞ ) ⎛ ∂x ⎞ ) 2 ⎛ ∂x ⎞ ) 2 m = ⎜ ⎟ mr2 + ⎜ ⎟ mu2 + ⎜ ⎟ m Ω + ⎜ ⎟ mi ⎝ ∂r ⎠ ⎝ ∂u ⎠ ⎝ ∂Ω ⎠ ⎝ ∂i ⎠ 2 x
Podívejme se, jakým dílem, z celkové střední chyby souřadnice x,
přispěje,
při
zanedbání
středních
chyb
všech
ostatních
veličin, střední chyba argumentu deklinace. Zaveďme několik předpokladů.
Hodnoty
rektascenze
výstupního
uzlu
dráhy
a
sklonu dráhy družice můžeme pro účely našeho odhadu brát jako konstanty.
Jedinou
proměnnou
veličinou
je
pro
nás
argument
deklinace. V závislosti na jeho hodnotě závorka v parciální derivaci dosahuje v druhé mocnině hodnoty až 0,908. Dále pro účely
tohoto
odhadu
postačí
vzít
jako
velikost
průvodiče
hodnotu r = 26560 km. Odhad střední chyby argumentu deklinace je mu = 0,00748 rad. A při maximální velkosti závorky ve výrazu parciální derivace je hodnota, kterou přispívá střední chyba argumentu
deklinace
k celkové
velikosti
střední
chyby
souřadnice x, rovna 189,3 km. Vlivy středních chyb ostatních veličin jsou i o několik řádů menší.Shrnu-li tyto výsledky, musím konstatovat, že střední chyby souřadnic družic určené jak na základě interpolace v terestrických souřadnicích tak i na základě interpolace v Keplerových dráhových elementech jsou velmi velké.
47
Pro porovnání původních dat s interpolovanými jsem použil jako
demonstračního
16),
ale
podobně
příkladu by
argumentu
vypadaly
i
deklinace
grafy
(obr.
srovnání
č.:
ostatních
Keplerových dráhových elementů a terestrických souřadnic.
Argument deklinace - u původní data
interpolované hodnoty
180 150 120
[o]
90 60 30 0
60
60
60
30
30
30
0
0
0 385
380
375
370
365
360
355
350
345
340
335
50
45
40
35
30
25
20
15
10
5
0
660
90
655
90
650
90
645
120
640
120
635
120
630
150
625
150
620
150
615
180
610
180
600
500
400
300
200
100
0
180
čas [počátek 50537.0625 MJD=1, 1 dílek po 15 minutách] obr. č.: 16
Pro výpočet vyrovnaných koeficientů interpolačních funkcí jsem použil statistický soubor dat z celého týdne, tedy všech 660 časových okamžiků (proč ne 672 je vysvětleno v kapitole [1.6]).
Interpolované
vyrovnaných
hodnoty
koeficientů
spočtené
interpolačních
na
základě
funkcí
mají
takto tu
vlastnost, že se na okrajích takového statistického souboru velmi
odchylují
s původními
od
jejich
hodnotami
velmi
původních dobře
hodnot.
Naproti
korespondují
tomu
uprostřed
48
tohoto statistického souboru (viz. obr. č.: 16). To znamená, že pro výpočet koeficientů interpolačních funkcí není vhodné použít všechna data z celého GPS týdne najednou. Mnohem přesnějších výsledků bychom pravděpodobně dosáhli použitím menšího statistického souboru dat, tzv. výběrového statistického souboru dat, který by zahrnoval kratší časové období. Tento výběrový soubor dat by obsahoval lichý počet členů (určení jejich počtu by záviselo na dalších výpočtech) a koeficienty interpolačních funkcí spočtené na základě takového souboru bychom přiřadili k prostřednímu bodu tohoto výběrového souboru.
Výpočet
bychom
například
provedli
postupně
pro
všechna data jednoho GPS týdne a získali bychom tabelované hodnoty
koeficientů
interpolačních
funkcí
pro
každý
časový
okamžik daného GPS týdne. Úloha by tedy přešla na interpolaci v koeficientech interpolačních funkcí. V těchto koeficientech by
bylo
patrně
snazší
interpolovat
a
následně
je
pak
extrapolovat v čase a na základě extrapolovaných koeficientů bychom
stanovili
mnohem
přesnější
extrapolované
terestrické
souřadnice, respektive Keplerovy dráhové elementy. Ještě
lepších
výsledků
bychom
patrně
dosáhli
použitím
metody kolokace nejmenších čtverců. Jedná se o metodu, kdy hledáme
korelační
vztahy
mezi
jednotlivými
hodnotami
statistického souboru a na jejich základě pak stanovíme jistou váhovou matici, kterou zapojíme do výpočtů vyrovnání. Protože se tato metoda zabývá pouze signálem, bylo by nejprve nutné odfiltrovat pozadí časových závislostí, tj. trendy. Potom by tato metoda měla dávat ještě lepší výsledky než postup, který jsem naznačil v předchozím odstavci.
49
[4] POUŽITÁ LITERATURA [1]
Kosmická geodézie a kosmická geodynamika (Milan Burša, Jan Kostelecký)
[2]
Geodetická astronomie 10 (Josef Kabeláč, Jan Kostelecký)
[3]
Dynamika umělých družic v tíhovém poli Země (Milan Burša, Georgij Karský, Jan Kostelecký)
[4]
Teorie chyb a vyrovnávací počet (Miroslav Hampacher, Vladimír Radouch)
[5]
Vyšší geodézie (Souřadné soustavy) (Miloš Cimbálník)
[6]
Vyšší geodézie 2 (Leoš Mervart, Miloš Cimbálník)
[7]
Numerical recipies in C (Teukolsky and al.)
[8]
Základy GPS a jeho praktické aplikace (Otakar Švábenský, Jan Fixel, Josef Weigel)
50
[5] PŘÍLOHY [5.1]
Příloha č.: 1 – Argumenty pro řešení nutace
ARG1 0., 0., 0., 0., -2., 0., 0., 0., 0., 1., 0., 1., 0.,-1., 0., 0., 2., 0., 0., 2., 0., 2., 0., 0., 1., 0., 0., 0., 1., 0., 1., 0., -1., 0., 0., 0., 1., 0., -1., 0., -1., 0., 1., 0.,
0., 0., 0., 0., 2., 0., 2.,-2., 0., 0., 2.,-2., 2.,-2., 2.,-2., 0.,-2., 0., 0., 2.,-2., 2., 0., 0., 0., 2., 0., 2., 0., 0.,-2., 2., 0., 0., 2., 0., 0., 0., 0., 2., 2., 2., 0.,
1., -171996., -174.2, 92025., 8.9, 2., 2062. , .2 , -895. , .5 , 1., 46. , .0 , -24. , .0 , 2., -13187. ,-1.6 , 5736. ,-3.1, 0., 1426. ,-3.4 , 54. ,-0.1, 2., -517. , 1.2 , 224. ,-0.6, 2., 217. ,-0.5 , -95. , 0.3, 1., 129. , 0.1 , -70. , 0.0, 0., 48. , 0.0 , 1. , 0.0, 0., 17. ,-0.1 , 0. , 0.0, 2., -16. , 0.1 , 7. , 0.0, 2., -2274. ,-0.2 , 977. ,-0.5, 0., 712. , 0.1 , -7. , 0.0, 1., -386. ,-0.4 , 200. , 0.0, 2., -301. , 0.0 , 129. ,-0.1, 0., -158. , 0.0 , -1. , 0.0, 2., 123. , 0.0 , -53. , 0.0, 0., 63. , 0.0 , -2. , 0.0, 1., 63. , 0.1 , -33. , 0.0, 1., -58. ,-0.1 , 32. , 0.0, 2., -59. , 0.0 , 26. , 0.0, 1., -51. , 0.0 , 27. , 0.0
ARG2 2., 0.,-2., 0., -2., 0., 2., 0., 1.,-1., 0.,-1., 0.,-2., 2.,-2., 2., 0.,-2., 0., 0., 0., 2.,-2., 0., 1., 0., 0., 0.,-1., 0., 0., -2., 0., 0., 2., 0.,-1., 2.,-2., 2., 0., 0.,-2., 0., 1., 2.,-2., 1., 0., 0.,-1., 2., 1., 0.,-2., 0., 0.,-2., 2., 0., 1.,-2., 2., 0., 1., 0., 0., -1., 0., 0., 1., 0., 1., 2.,-2., 0., 0., 2., 2., 2., 0., 0., 0., 1., 0., 2.,-2., 2., 0., 2., 0., 0., 0., 2., 0., -1., 0., 2., 0., -1., 0., 0., 2., 1., 0., 0.,-2., -1., 0., 2., 2., 1., 1., 0.,-2., 0., 1., 2., 0., 0.,-1., 2., 0., 1., 0., 2., 2., 1., 0., 0., 2., 2., 0., 2.,-2., 0., 0., 0., 2., 0., 0., 2., 2., 1., 0., 2.,-2., 0., 0., 0.,-2., 1.,-1., 0., 0., 2., 0., 2., 0., 0., 1., 0.,-2., 1., 0.,-2., 0.,
0., 2., 0., 1., 1., 0., 1., 1., 1., 1., 1., 1., 0., 0., 1., 0., 2., 1., 0., 2., 0., 2., 2., 0., 1., 1., 1., 1., 0., 2., 2., 2., 0., 2., 1., 1., 1., 1., 0., 1., 0., 0.,
11. -3. -3. -2. 1. -22. -15. -12. -6. -5. 4. 4. -4. 1. 1. -1. 1. 1. -1. -38. 29. 29. -31. 26. 21. 16. -13. -10. -7. 7. -7. -8. 6. 6. -6. -7. 6. -5. 5. -5. -4. 4.
, 0. , , 1. , , 0. , , 1. , , 0. , , 0. , , 9. , , 6. , , 3. , , 3. , , -2. , , -2. , , 0. , , 0. , , 0. , , 0. , , 0. , , 0. , , 0. , , 16. , , -1. , , -12. , , 13. , , -1. , , -10. , , -8. , , 7. , , 5. , , 0. , , -3. , , 3. , , 3. , , 0. , , -3. , , 3. , , 3. , , -3. , , 3. , , 0. , , 3. , , 0. , , 0. ,
0., 0., 0., 1., 1., 1., 0., 0., 1., 0., 2., 0., 1.,-1., 2., 0., -1.,-1., 2., 2., -2., 0., 0., 0., 3., 0., 2., 0., 0.,-1., 2., 2., 1., 1., 2., 0., -1., 0., 2.,-2., 2., 0., 0., 0., 1., 0., 0., 0., 3., 0., 0., 0., 0., 0., 2., 1., -1., 0., 0., 0., 1., 0., 0.,-4., -2., 0., 2., 2., -1., 0., 2., 4., 2., 0., 0.,-4., 1., 1., 2.,-2., 1., 0., 2., 2., -2., 0., 2., 4., -1., 0., 4., 0., 1.,-1., 0.,-2., 2., 0., 2.,-2., 2., 0., 2., 2., 1., 0., 0., 2., 0., 0., 4.,-2., 3., 0., 2.,-2., 1., 0., 2.,-2., 0., 1., 2., 0., -1.,-1., 0., 2., 0., 0.,-2., 0., 0., 0., 2.,-1., 0., 1., 0., 2., 1., 0.,-2.,-2., 0.,-1., 2., 0., 1., 1., 0.,-2., 1., 0.,-2., 2., 2., 0., 0., 2., 0., 0., 2., 4., 0., 1., 0., 1.,
0., 0., 0., 2., 2., 1., 2., 2., 2., 1., 1., 2., 0., 2., 2., 0., 2., 2., 0., 2., 1., 2., 2., 0., 1., 2., 1., 2., 2., 0., 1., 1., 1., 2., 0., 0., 1., 1., 0., 0., 2., 0.,
-4. -3. 3. -3. -3. -2. -3. -3. 2. -2. 2. -2. 2. 2. 1. -1. 1. -2. -1. 1. -1. -1. 1. 1. 1. -1. -1. 1. 1. -1. 1. 1. -1. -1. -1. -1. -1. -1. -1. 1. -1. 1.
, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
0. 0. 0. 1. 1. 1. 1. 1. -1. 1. -1. 1. 0. -1. -1. 0. -1. 1. 0. -1. 1. 1. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,
51
[5.2]
Příloha č.: 2 –
Seznam řídících a
monitorovacích stanic hlavní řídící stanice:
Colorado Springs (USA)
pozemní řídící stanice: Ascension Island (jižní Atlantik) Diego Garcia (Čagoské ostrovy, Indický oceán) Kwajalwin (Marshallovo souostroví, Tichý oceán) monitorovací stanice:
Hawai (Tichý oceán) Ascension Island (jižní Atlantik) Diego Garcia (Čagoské ostrovy, Indický oceán) Kwajalein (Marshallovo souostroví, Tichý oceán) Colorado Springs (USA)
52
[5.3]
Příloha č.: 3 – Programy InterKPL a InterXYZ
V této
příloze
popíši
oba
programy,
které
jsem
pro
výpočetní zajištění této práce napsal. Kapitola je rozdělena do dvou podkapitol. Každá z nich je věnována jednomu ze dvou programů. Nyní bych chtěl napsat několik řádek obecně o programech jako o jednom celku. Oba jsou psány v jazyce C++. Jak již bylo výše
v práci
řečeno,
při
psaní
programů
jsem
použil
interaktivní objektově orientované prostředí C++ Builder 3.0. Oba programy jsou psány jako aplikace pod systémy WIN95, WIN98 a WIN NT. Vstupními formáty souborů jsou soubory SP3 pro vstup terestrických souřadnic družic a ERP pro vstup souřadnic pólu a korekce (UT1-UTC). Některé výpočetní bloky byly převzaty ze zdrojových textů publikovaných
Astronomickým
(Astronomical
institute,
institutem
University
univerzity
of
Bern).
v Bernu
Nicméně
jsem
tyto bloky musel převést z jazyka FORTRAN do jazyka C++. To obnášelo
alespoň
programovacího
základní
jazyka.
Na
znalosti
tyto
bloky
syntaxe
bude
tohoto
upozorněno
dále
v textu. Celkový
objem
zdrojových
kódů
obnáší
4726
programových
řádek. Jsou to poměrně rozsáhlé programy a proto není možné v této
práci
jako
přílohu
uvést
jejich
vytištěné
zdrojové
texty. Tyto texty budou na přiloženém CD, které bude obsahovat i ostatní části této diplomní práce v digitální podobě. Obsah přiloženého
CD
je
uveden
v příloze
č.:4.
Stejně
tak
nemá
příliš opodstatnění uvést popis jednotlivých funkcí, metod a procedur, protože jich je mnoho a pak také jsem programy psal tak,
že
patrné
je
z názvů
k jakému
jednotlivých
účelu
slouží.
procedur, V této
metod
příloze
a
budou
funkcí tedy
uvedeny jen některé zajímavé funkce a metody, které nezaberou tolik místa. A nyní již k jednotlivým programům.
53
[5.3.1]
Program InterKPL
Tento program řeší výpočet Keplerových dráhových elementů a jejich následnou interpolaci dle vzorců uvedených v kapitole [1.9].
Po
spuštění
se
otevře
vlastní
okno
programu
(obr.
č.:17). U horního okraje je menu s položkami, jejichž funkci vysvětlím níže. Pod tímto menu je komunikační okno, pomocí něhož pak program při výpočtu signalizuje jeho stav. Pod tímto
obr. č.: 17
oknem je další okno, které má název „Asociované soubory:“. Do tohoto
okna
se
vypisují
soubory,
které
byly
určeny
jako
vstupní. Pod tímto oknem jsou tři stavové řádky. První ukazuje stav právě probíhajících výpočtů. Druhý ukazuje číslo družice pro kterou se právě výpočet provádí a dále délku ukončeného výpočtu. Poslední stavová řádka ukazuje aktuální adresář. Před samotným výpočtem je třeba provést několik nastavení. Za
prvé
výpočet
je
nutné
provést.
nastavit To
se
číslo
provádí
družice, pomocí
pro menu
kterou
chci
[Program
Æ
Nastavení Æ Číslo družice]. Implicitně je nastavena družice PRN1. Dále se musí zadat soubor se souřadnicemi pólu a korekcí (UT1-UTC). To se provede v menu [Program Æ Nastavení Æ Soubor
54
*.erp].
Dále
je
třeba
nastavit,
kolik
dní
bude
vzato
do
výpočtu. Implicitně nastaveno je 7, tedy celý týden. Klávesami F2-F7
nebo
přímo
v položce
Æ
[Soubory
Min.
počet
asoc.
souborů] je možné tento údaj měnit. Dalším nutným krokem je tyto soubory asociovat. Je možné tento úkon provést kombinací kláves
CTRL+A
a
nebo
v
menu
Æ Asociovat]. Každý
[Soubory
soubor se musí asociovat samostatně. Pokud chceme odstranit nějaký
soubor
z našeho
výběru,
můžeme
jej
odstranit
dvojím
kliknutím myši na název souboru s cestou v panelu „Asociované soubory:“ a nebo v menu [Soubory Æ Deasociovat]. Posledním úkonem,
který
před
vlastním
výpočtem
je
třeba
provést,
je
nastavení přibližných hodnot amplitud, fázových posunů, period a trendů. To se provádí pomocí menu [Program Æ Nastavení Æ Přibližné
koeficienty]
a
nebo
dvojím
kliknutím
myší
na
prostřední stavový řádek dole. Implicitně jsou zde nastaveny přibližné hodnoty těchto veličin pro družici PRN1 a GPS týden 899. Kliknutím na volnou plochu programu se zadávání ukončí. Nyní
můžeme
spustit
samotný
výpočet.
Provedeme
kontrolu
asociovaných souborů pomocí menu [Soubory Æ Kontrola vstupní konfigurace]
nebo
výpočet
v menu
CTRL+V.
Výsledný
kombinací
[Výpočet
Æ
výpočetní
kláves
CTRL+K
Spustit]
nebo
protokol
lze
a
pak
spustíme
kombinací pak
uložit
kláves pomocí
kláves CTRL+U nebo v menu [Soubory Æ Záznamy výpočtů Æ Uložit záznam výpočtů]. Ještě je dobré říci, že program je vybaven nápovědou v podobě tzv. HINTů, tj. pokud kurzor myši zůstane na části programu déle, objeví se vedle něj okénko s průvodním textem.
A
poslední
výhodné
v menu
poznámka
[Soubory
Æ
k tomuto Změna
programu
adresáře]
je,
nebo
že
je
klávesami
CTRL+Z nastavit adresář vstupních souborů *.sp3. Nyní bych rád uvedl příklad převzatého výpočetního bloku, který jsem přepsal z jazyka FORTRAN do jazyka C++. Jedná se výpočet nutační matice. Tedy původní zdrojový text v jazyce FORTRAN má tento formát:
55 C* SUBROUTINE NUTN20(XTDB,NUT) CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC CC C*
NAME
:
NUTN20
CALL NUTN20(XTDB,NUT) PURPOSE
:
COMPUTATION OF NUTATION-MATRIX FOR EPOCH XTDB (JUL. DATE IN BARYCENTRIC DYNAMICAL TIME) SEE "USNO CIRCULAR NO 163",1981 (IAU RESOLUTIONS) EDITED BY KAPLAN
-> OLD VARIABLE XMOD IS POSSIBLE TOO, BECAUSE DIFF. < 1 MIN CAUSES DIFFERENCIES < 0.1 MAS PARAMETERS : IN : OUT :
XTDB NUT
SR CALLED
:
GETCPO, LITPOL
REMARKS
:
---
AUTHOR
:
E. BROCKMANN
VERSION
:
3.3
CREATED
:
92/05/05 09:20
CHANGES
:
25-AUG-92 : USE CORRECT SIN AND COS FUNCTIONS TO COMPUTE ROTATION MATRIX 20-OCT-92 : USE CORRECT DEPS AND DPSI INSTEAD OF DMUE AND DNUE. 29-OCT-93 : SF: UPDATE HEADER INFORMATION 19-APR-94 : RW: CPO-MODEL INCLUDED
COPYRIGHT 1992
:
ASTRONOMICAL INSTITUTE UNIVERSITY OF BERNE SWITZERLAND
: EPOCH IN BARYCENTRIC DYNAMICAL TIME : NUTATION - MATRIX
LAST MODIFIED :
R*8 R*8(3,3)
19-APR-94
IMPLICIT REAL*8 (A-H,O-Z) REAL*8 NUT(3,3),ARG1(9,22),ARG2(7,84),ALP(5) REAL*8 TT(2),TTEPS(2),TTPSI(2) C C DATA C ---DATA PI/3.141592653589793D0/,R/1296000.D0/ C C 22 COEFFICIENTS FOR THE STONGEST NUTATIONPERIODS (TIMEDEPENDENT) C C UNIT: ARCSEC/10000 C C ROH = PI/648000.D0 C C TIME INTERVAL (IN JUL. CENTURIES) BETWEEN XTDB AND J2000.0 TU =(XTDB-51544.5)/36525.D0 C TU4 = MOD JUL. DATUM OF REFERENCE EPOCH ( XMOD = XTDB FOR EPS) TU4=TU C C FUNDAMENTAL ARGUMENTS (IN RAD) ALP(1)=(485866.733+(1325.D0*R+715922.633)*TU+31.31*TU*TU+ 1 .064*TU*TU*TU)*ROH ALP(2)=(1287099.804+(99.D0*R+1292581.224)*TU-0.577*TU*TU1 .012*TU*TU*TU)*ROH ALP(3)=(335778.877+(1342.D0*R+295263.137)*TU-13.257*TU*TU+ 1 .011*TU*TU*TU)*ROH ALP(4)=(1072261.307+(1236.D0*R+1105601.328)*TU-6.891*TU*TU+ 1 .019*TU*TU*TU)*ROH ALP(5)=(450160.28-(5.D0*R+482890.539)*TU+7.455*TU*TU+ 1 .008*TU*TU*TU)*ROH C C EPS-ZERO EPOCH J2000 (RAD) EPS=84381.448-46.815*TU4-0.00059*TU4*TU4+0.001813*TU4*TU4*TU4 EPS=EPS*ROH C
56 DPSI=0.D0 DEPS=0.D0 DO 10 J=1,22 PHI = 0.D0 DO 5 I=1,5 PHI = PHI + ARG1(I,J)*ALP(I) 5 CONTINUE PHI = DMOD (PHI,2.D0*PI) C NUTATION IN ARCSEC*10000 DPSI = DPSI+(ARG1(6,J)+ARG1(7,J)*TU)*DSIN(PHI) DEPS = DEPS+(ARG1(8,J)+ARG1(9,J)*TU)*DCOS(PHI) 10 CONTINUE DO 20 J=1,84 PHI = 0.D0 DO 15 I=1,5 PHI = PHI + ARG2(I,J)*ALP(I) 15 CONTINUE PHI = DMOD (PHI,2.D0*PI) C SERIES FOR NUTATION IN ARCSEC*10000 IN LONGITUDE AND OBLIQUITY DPSI = DPSI+ARG2(6,J)*DSIN(PHI) DEPS = DEPS+ARG2(7,J)*DCOS(PHI) 20 CONTINUE C C DPSI AND DEPS IN RADIANS DEPS=DEPS*ROH/10000.D0 DPSI=DPSI*ROH/10000.D0 C C CELESTIAL POLE OFFSETS IN RADIANS CALL GETCPO(XTDB,TT,TTEPS,TTPSI) CALL LITPOL(2,1,TT,TTEPS,XTDB,DDEPS) CALL LITPOL(2,1,TT,TTPSI,XTDB,DDPSI) DEPS=DEPS+DDEPS DPSI=DPSI+DDPSI C C TRUE OBLIQUITY EPSTRU EPSTRU=EPS+DEPS C C SINE AND COSINE SINDP=DSIN(DPSI) COSDP=DCOS(DPSI) SINEM=DSIN(EPS) COSEM=DCOS(EPS) SINET=DSIN(EPSTRU) COSET=DCOS(EPSTRU) C C ROTATION MATRIX NUT(1,1)= COSDP NUT(2,1)= COSET*SINDP NUT(3,1)= SINET*SINDP NUT(1,2)=-COSEM*SINDP NUT(2,2)= COSEM*COSET*COSDP + SINEM*SINET NUT(3,2)= COSEM*SINET*COSDP - SINEM*COSET NUT(1,3)=-SINEM*SINDP NUT(2,3)= SINEM*COSET*COSDP - COSEM*SINET NUT(3,3)= SINEM*SINET*COSDP + COSEM*COSET C RETURN END
Za
řádkem
C
DATA
následuje
seznam
argumentů
pro
výpočet
lineární kombinace Delauneyových proměnných. Tento seznam jsem z ukázky odebral, neboť osahuje jen samá čísla a byl několik stánek
dlouhý.
Seznam
těchto
argumentů
je
uveden
v příloze
č.:1. A
mnou
napsaný
zdrojový
text,
který
řeší
stejnou
problematiku, tedy počítá nutační matici, je následující:
57
Mat<>MaticeNutace(double TGPS) { double R=1296000; double roh,eps,epstru; double sindp,cosdp,sinem,cosem,sinet,coset; double tu,tu4; //cas mezi J2000.0 a aktual. epochou ve stoletich double* alp; alp = new double[5]; double dpsi, deps; double phi; Mat<>Nutace; Nutace.reset(3,3); double TDT=TGPS+(19+32.184)/86400; tu=(TDT-51544.5)/36525; tu4=tu; roh=M_PI/648000; alp[0]=(485866.733+(1325*R+715922.633)*tu+31.31*tu*tu+0.064*tu*tu*tu)*roh; alp[1]=(1287099.804+(99*R+1292581.224)*tu-0.577*tu*tu-0.012*tu*tu*tu)*roh; alp[2]=(335778.877+(1342*R+295263.137)*tu-13.257*tu*tu+0.011*tu*tu*tu)*roh; alp[3]=(1072261.307+(1236*R+1105601.328)*tu-6.891*tu*tu+0.019*tu*tu*tu)*roh; alp[4]=(450160.28-(5*R+482890.539)*tu+7.455*tu*tu+0.008*tu*tu*tu)*roh; // epsilon 0 v radianech eps=(84381.448-46.815*tu4-0.00059*tu4*tu4+0.001813*tu4*tu4*tu4)*roh; dpsi=0; deps=0; for(int j=1;j<=22;j++){ phi=0; for(int i=1;i<=5;i++){ phi=phi+ARG1(i,j)*alp[i-1]; } phi=phi-INT(phi/(2*M_PI))*(2*M_PI);// zbytek po deleni //nutace v arcsec*10000 dpsi=dpsi+(ARG1(6,j)+ARG1(7,j)*tu)*sin(phi); deps=deps+(ARG1(8,j)+ARG1(9,j)*tu)*cos(phi); } for(int j=1;j<=84;j++){ phi=0; for(int i=1;i<=5;i++){ phi=phi+ARG2(i,j)*alp[i-1]; } phi=phi-INT(phi/(2*M_PI))*(2*M_PI);// zbytek po deleni //nutace v arcsec*10000 v zemepisne delce a "sikmosti"} dpsi=dpsi+ARG2(6,j)*sin(phi); deps=deps+ARG2(7,j)*cos(phi); } //dpsi a deps v radianech dpsi=dpsi*roh/10000; deps=deps*roh/10000; //celkova nutace epstru=eps+deps; //siny a cosiny pro matici rotace sindp=sin(dpsi); cosdp=cos(dpsi); sinem=sin(eps); cosem=cos(eps); sinet=sin(epstru); coset=cos(epstru); //nutacni matice rotace Nutace(1,1)=cosdp; Nutace(2,1)=coset*sindp; Nutace(3,1)=sinet*sindp; Nutace(1,2)=-cosem*sindp; Nutace(2,2)=cosem*coset*cosdp+sinem*sinet; Nutace(3,2)=cosem*sinet*cosdp-sinem*coset; Nutace(1,3)=-sinem*sindp; Nutace(2,3)=sinem*coset*cosdp-cosem*sinet; Nutace(3,3)=sinem*sinet*cosdp+cosem*coset; delete[]alp; return Nutace; }
58
[5.3.2]
program InterXYZ
Program InterXYZ řeší po výpočetní stránce problematiku interpolace v terestrických souřadnicích družic. Jeho ovládání je podobné jako u programu InterKPL, ale je jednodušší. Po spuštění se objeví opět samostatné okno tohoto programu. Toto okno je až na název totožné s oknem programu InterKPL. Asociace, změna adresáře, ukládání protokolů, deasociace souborů a nastavení počtu dnů vzatých do výpočtů je naprosto stejné
jako
u
programu
InterKPL.
Změnilo
se
jen
zadávání
přibližných koeficientů a nastavení PRN počítané družice. Dále odpadlo nastavení souboru se souřadnicemi pólu a korekcí (UT1UTC). Při výpočtech postupujeme stejně jako v kapitole [5.3.1] až na malé rozdíly, které jsou: zadání přibližných koeficientů interpolačních [Soubory pomocí IterKPL,
Æ
funkcí
Přibližné
klávesové jsou
i
(viz.
koeficienty],
zkratky zde
kapitola
CTRL+P.
přednastaveny
[1.8]) nebo
Podobně
je
je
možné
jako
přibližné
zde u
v menu vyvolat
programu
hodnoty
těchto
koeficientů pro družici PRN1 a GPS týden 899. Zadání čísla družice
je
v menu
[Soubory
klávesové zkratky CTRL+C.
Æ
Číslo
družice],
nebo
pomocí
59
[5.4]
Příloha č.: 4 – Obsah přiloženého CD
Adresář
\Data
:
obsahuje data pro GPS týdny 899 a 900. Tj. terestrické družic
a
korekce
souřadnice
souřadnice
(UT1-UTC)
pólu
a
to
a jak
finální tak i predikované. \Programy
:
obsahuje
spustitelné
programy
InterKPL a InterXYZ, zajišťující po výpočetní stránce problematiku této diplomové práce. \Různé
:
obsahuje
různé
soubory TXT,
ve formátech EXCEL 97 a
vzniklé
tématu.
doprovodné
během
Tyto
jakýchkoliv
zpracovávání
soubory
jsou
grafických
bez
úprav
a
obsahují hrubé dokumenty. \Texty
:
obsahuje text diplomové práce ve formátu WORD 97 a RTF.
\Zdrojové kódy :
obsahuje
všechny
zdrojové texty
programů, které jsem programoval v jazyce
C++,
některých
zdrojové
výpočetních
v jazyce
FORTRAN
zdrojové
texty
souborů
gmatvec
s maticemi, Ing.
Aleš
texty
které Čepek,
bloků
a
konečně
hlavičkových pro napsal CSc
práci doc.
z katedry
Mapování a kartografie a které mi poskytl.
60
[5.5]
Příloha č.: 5 – predikovaných
Srovnání oficiálních souřadnic
družice
s jejich finálními hodnotami
GPS
PRN1
61
62
63
64
