VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ FAKULTA STAVEBNÍ INFORMATIKA II MODUL 03 ALGORITMIZACE A PROGRAMOVÁNÍ V INŽENÝRSKÝCH ÚLOHÁCH

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ FAKULTA STAVEBNÍ

INFORMATIKA II MODUL 03 ALGORITMIZACE A PROGRAMOVÁNÍ V INŽENÝRSKÝCH ÚLOHÁCH

STUDIJNÍ OPORY PRO STUDIJNÍ PROGRAMY S KOMBINOVANOU FORMOU STUDIA

Informatika II • Modul 03 • Algoritmizace a programování

© Jiří Macur, Miroslav Menšík

2

Obsah

OBSAH OBSAH

3

VSTUPNÍ INFORMACE K MODULU

5

Cíle

5

Požadované znalosti

5

Doba potřebná ke studiu

5

Klíčová slova

5

Použitá terminologie

5

Metodický návod na práci s textem

5

1

PRŮBĚH FUNKCE

6

1.1

List s parametry úlohy a výstupním grafem

6

1.2

Implementace

7

1.3

Pokročilejší techniky

8

2

NUMERICKÁ INTEGRACE A DERIVACE

11

3

POČÁTEČNÍ ÚLOHA

16

3.1

Implementace výpočtu trajektorie

17

3.2


18

3.3

Interaktivní prvek posuvník

20

3.4

Použití pole v počáteční úloze

22

3.5

Balistická křivka

24

4

SIMULACE POHYBU PRUŽNÉHO LANA

4.1

Fyzikální model

26

4.2

Tlumení

27

4.3

Implementace

28

4.4

Komentář

30

4.5

Zastavení simulačního programu

32

3

26

Informatika • Modul 03 • Algoritmizace a programování 4.6

Úkoly

34

5

ZPRACOVÁNÍ NMEA DAT Z MĚŘICÍHO PŘÍSTROJE

5.1

Satelitní čas

37

5.2

Transformace souřadnic

37

5.3

Implementace zpracování dat

40

5.4

Komentář

40

4

35

Řešené příklady

Vstupní informace k modulu Cíle Modul obsahuje vybrané řešené příklady inženýrských úloh v prostředí jazyka Visual Basic a objektových knihoven vybavení MS Office. Cílem textu je ukázat možnosti jazyka a objektových aplikací při řešení a vizualizaci standardních problémů.

Požadované znalosti K přečtení a porozumění tomuto textu jsou potřebné znalosti uvedené v modulu CU01_M01 a CU01_M02 a práce v prostředí VBA.

Doba potřebná ke studiu Je velmi individuální podle znalostí čtenáře. Při pozorném studiu programů, resp. jejich rozvíjení o další funkcionalitu, se může jednat o cca 15 hod.

Klíčová slova Algoritmizace, programovací jazyk Visual Basic, strukturované programování, objektově orientované programování, objektový model aplikace, knihovny objektů, implementace inženýrských úloh.

Použitá terminologie Diskutovaná problematika se rychle vyvíjí, proto je česká terminologie značně neustálená. V textu je většinou uváděna i ekvivalentní anglická terminologie spolu s příslušnými zkratkami.

Metodický návod na práci s textem Programy by měly být čtenářem do hloubky pochopeny. Jsou krátké a reprezentují základní postupy při využívání znalostí z předcházejících modulů. Rozhodně není vhodné přeskakovat jednotlivé příkazy a části programů s tím, že budou pochopeny později. Pokud čtenáři není zřejmá uvedená algoritmická konstrukce nebo syntaxe příkazu, měl by si vyhledat potřebné informace v bohatých odkazech internetu nebo požádat o konzultaci vyučujícího.

5

Informatika • Modul 03 • Algoritmizace a programování

1

Průběh funkce

MS Excel poskytuje dostatek nástrojů pro zobrazení průběhu funkce. Standardně lze do jednoho sloupce vložit řadu pro nezávisle proměnnou, do buňky jiného sloupce můžeme vložit formuli obsahující funkci s odkazem na buňku nezávisle proměnné. Formuli zkopírujeme do potřebné oblasti a získáme řadu závisle proměnné. Nad oběma řadami pak lze vytvořit bodový graf X–Y. Postup však můžeme také automatizovat jednoduchým programem, na němž lze demonstrovat několik používaných technik.

1.1

List s parametry úlohy a výstupním grafem

Návrh listu vytvoříme interaktivně například takto:

Obr. 1.1 – Návrh interaktivního listu pro automatický průběh funkce List nazveme "funkce" a druhý list v sešitu nazveme "data". Ostatní listy v sešitu odstraníme. Do prvního listu vložíme bodový graf, v interaktivním formuláři vložíme libovolné buňky (alespoň dvě) tak, aby se vytvořila v grafu datová řada:

6


1.2

Implementace

Nyní vytvoříme v programovém modulu jednoduchý program:

Obr. 1.2 – Implementace průběhu funkce Program má jasně definované sekce: 7


  

Načtení dat z listu do proměnných programu Výpočet pomocí cyklu hodnot nezávisle a závisle proměnné uživatelské funkce a vložení vypočtených hodnot do patřičných buněk v listu "data". Nastavení vlastností grafu – titulek (podobjekt ChartTitle) a rozsahy buněk tvořících datovou řadu tak, jak je popsáno v modulu CU01_M02.

Nakonec do interaktivního listu vložíme tlačítko pro spuštění programu (karta Vývojář – Vložit – Prvky formuláře) a přiřadíme mu náš podprogram Funkce. V listu "data" se budou při běhu programu buňky přepisovat, v grafu se automaticky použije nově vypočtená oblast.

1.3


Naše malá aplikace má jednu nevýhodu. Uživatel musí při změně zobrazované funkce opravit rovněž funkci f() v programovém modulu. Komfortnější by byla možnost definovat funkci pouze v buňce B4 – všechno ostatní by bylo provedeno programem. Tento požadavek není ve většině programovacích jazyků snadno splnitelný – ve skutečnosti se jedná o provádění příkazů, které v době překladu ještě nejsou známy. Jsou vytvořeny až za běhu programu. V objektu Application však máme k dispozici speciální metodu Evaluate, která pracuje přesně podle výše uvedeného požadavku: vyhodnotí výraz, který není součástí žádného příkazu – je uložen pouze v podobě řetězce, který předáme metodě v podobě parametru. Obsah řetězce musí být z hlediska syntaxe jazyka VB samozřejmě správný, jinak dojde k chybě za běhu programu. Naši funkci f tedy můžeme modifikovat takto: Function f(ByVal x As Single) As Single Dim strfunkce As String strfunkce = Worksheets("funkce").Cells(4, 2).Value strfunkce = Replace(strfunkce, "x", CStr(x)) f = Application.Evaluate(strfunkce) End Function

Nejprve vytvoříme řetězec, který obsahuje výraz uložený v buňce B4. V tomto výrazu musíme zaměnit symbolický argument "x" zadané funkce za hodnotu nezávisle proměnné, z níž chceme vypočítat odpovídající funkční hodnotu. K tomu nám poslouží řetězcová funkce Replace, která všechny výskyty písmena "x" v řetězci nahradí číslem (převedeným na řetězec) předaným do funkce v podobě hodnoty parametru. Metoda Evaluate pak vytvořený řetězec (v našem případě např. "sin(0.45678)/0.45678") vyhodnotí, jako by to byl výraz v programu. Poznámka Uvedený postup byste měli chápat spíše jako trik, který příliš často v praxi nepoužijeme. Je pomalý a v našem případě má i jiné nectnosti – např. v případě, že ve výrazu bude použito písmeno "x" i jinde, než ve významu argumentu (např. funkce exp(x)), funkce samozřejmě selže. Bude se totiž vyhodnocovat 8


nesmyslný výraz typu "e0.32p(0.32)". Tento problém můžeme vyřešit složitějším označením argumentu funkce (např. arg). Podstatnějším problémem je však národní interpretace desetinné čárky (resp. tečky). Podle konfigurace aplikace je třeba, aby metoda Evaluate vyhodnocovala čísla ve výrazu s desetinnou tečkou nebo čárkou; nevhodný znak vyvolá chybu. Pro uživatele našeho programu může být přítomnost dalšího listu "data" rušivá. Navíc, ukládání vypočtených hodnot do listu, odkud jsou opět načítány do grafu, může program zbytečně zpomalovat. Potřebné řady dat můžeme uložit do polí a ta předat grafu bez mezikroku. Musíme ovšem použít dynamická pole, neboť počet jejich prvků stanoví uživatel až v době běhu programu. Použili jsme typ cyklu For … Next, který je v uvedeném případě pohodlnější. Implementace programu se poněkud zjednoduší a program se urychlí:

Obr. 1.3 – Návrh s využitím pole a výpočtem funkce podle obsahu buňky

9


Nyní je program komfortní, zobrazuje průběh zadané funkce bez nutnosti dalších úprav. Překvapivá jednoduchost programu je daná efektivitou použitých knihoven, v tomto případě zejména třídou Chart. Poznámka Počet prvků v poli, které předáváme do datové řady grafu je omezen. Pokud bude bodů příliš mnoho, použije se pro tvorbu grafu pouze část z nich. Důvod nalezne čtenář v poznámkách o objektovém modelu grafu v modulu CU01_M02.

Obr. 1.4 – Jediný list aplikace poskytující zobrazení průběhu funkce

10


2 Numerická integrace a derivace Velké množství určitých integrálů nelze vyjádřit prostřednictvím elementárních funkcí, např. při integraci transcendentních funkcí. V takovém případě používáme k určení hodnoty integrálu přibližných metod, mezi něž patří tzv. numerická integrace. Při numerické integraci se snažíme nahradit integrál jiným druhem výpočtu, přičemž se snažíme zajistit, aby se získaná hodnota od skutečné hodnoty integrálu lišila co nejméně. Metod numerické integrace existuje velké množství. Výběr správné metody závisí na požadované přesnosti řešení, časové náročnosti a jiných parametrech. Mezi nejjednodušší metody numerické integrace patří např. obdélníková, lichoběžníková nebo Simpsonova metoda. Pro naše potřeby plně vyhoví jednoduchá lichoběžníková metoda:

Při lichoběžníkové metodě numerické integrace rozdělíme interval pomocí bodů xi, kde a = x0 < x1 < x2 < ... < xn = b na n stejně velkých podintervalů o velikosti h = (b – a)/n. Určitý integrál (plochu pod křivkou funkce) při použití této metody aproximujeme jako součet ploch lichoběžníků podle výše uvedeného obrázku:

Poznámka Výsledný vztah pro výpočet integrálu lze upravit na tvar , z dobrých důvodů však použijeme předcházející součet.

11


Pro výpočet integrálu použijeme program z obr. 1.2. Pro sečítání ploch lichoběžníků využijeme cyklus výpočtu bodů pro graf: For i = 1 To n poleX(i) = x poleY(i) = f(x) If i < 0 Then plocha = plocha + 0.5 * h * (f(x) + f(x + h)) x = x + h Next i ' -- výstup určitého integrálu Worksheets("funkce").Cells(8, 2).Value = plocha

Podmíněný příkaz v přičítání ploch zamezí přičtení plochy lichoběžníku, který leží již za horní mezí integrálu. Výstup programu pak může vypadat např. takto:

Obr. 2.1 – Pracovní list průběhu funkce s výpočtem určitého integrálu v zadaných mezích Všimněme si, že název argumentu jsme museli změnit kvůli zmiňovanému výskytu písmena "x" v názvu exponenciální funkce. Zadaná gaussova funkce má v mezích (-∞, ∞) hodnotu integrálu rovnu √π = 1.772453851. Vidíme, že funkce rychle konverguje k nule, protože již v mezích <-5, 5> je numericky vypočtený integrál prakticky totožný. Zvolený postup s výhodou využijeme pro zobrazení průběhu neurčitého integrálu, který můžeme chápat jako funkci jeho horní meze:

Dolní mez může být libovolná (samozřejmě funkce i integrál musí v uvedených mezích existovat) – její neurčitost je zahrnuta v integrační konstantě C. V programu podle obr. 1.2 provedeme následující úpravy: přidáme proměnnou typu pole, kterou deklarujeme stejně, jako pole bodů pro graf průběhu funkce.

12


Do cyklu přidáme příkazy s parciálními součty ploch lichoběžníků, které postupně uložíme do jednotlivých prvků nového pole např. takto: For i = 1 To n poleX(i) = x poleY(i) = f(x) poleI(i) = plocha plocha = plocha + 0.5 * h * (f(x) + f(x + h)) x = x + h Next i Proměnná poleI bude obsahovat funkci neurčitého integrálu zadané funkce. Do grafu přidáme interaktivně další datovou řadu:

Na rozsahu buněk pro řadu opět nezáleží, budou přepsány následujícími programovými řádky: Worksheets("funkce").ChartObjects(1).Chart.SeriesCollection(2).XValues = poleX Worksheets("funkce").ChartObjects(1).Chart.SeriesCollection(2).Values = poleI

Index v kolekci datové řady samozřejmě musí být pro novou řadu roven 2. Náš program nyní zobrazuje nejen průběh funkce, ale i průběh jejího neurčitého integrálu. Stejně jednoduše můžeme do grafu přidat další datovou řadu, která zobrazí průběh první derivace funkce. Vytvoříme další pole poleD a v cyklu programu ji naplníme hodnotami s aproximací:

Do grafu přidáme interaktivně další řadu a zobrazení vypočteného průběhu derivace funkce zajistíme programovými řádky: Worksheets("funkce").ChartObjects(1).Chart.SeriesCollection(3).XValues = poleX Worksheets("funkce").ChartObjects(1).Chart.SeriesCollection(3).Values = poleD

Všimněme si, že musíme přidat pro každou řadu rovněž pole hodnot nezávisle proměnné i v případě, že jsou tyto obsahy pro každou datovou řadu stejné. Implementace celého programu tedy bude následující:

13


Obr. 2.1 – Implementace zobrazení průběhu funkce, neurčitého integrálu a derivace

14


Pracovní list celé aplikace pak může vypadat takto:

V listu můžeme změnit nejen parametry zobrazované oblasti, ale i funkci samotnou. Pro příliš odlišné datové hodnoty derivace od hodnot funkce ("dramatický" průběh funkce) můžeme pro tuto datovou osu přesunout měřítko na vedlejší osu:

15


3

Počáteční úloha

V úlohách dynamiky obvykle hledáme vývoj systému v čase na základě pohybových rovnic. V jednoduchých případech hmotného bodu je představuje Newtonůw zákon síly: F = ma, kde F je vektor působící síly, m je hmotnost, a zrychlení. Tento pohled lze zobecnit zavedením vektoru stavu X, který jednoznačně a úplně určuje další vývoj systému. Pro hmotný bod je takovým vektorem stavu dvojice polohového vektoru a vektoru rychlosti. V rovinné úloze tedy bude mít stavový vektor čtyři položky X ≡ (x, y, vx, vy). Principem numerického řešení evoluce systému je: 1. Diskretizace času – stav určujeme pouze v posloupnosti časových okamžiků ti mezi časovými intervaly (kroky): Δti = ti+1 – ti. Velikost časového kroku se tradičně označuje písmenem h = Δti a je obvykle konstantní (intervaly jsou ekvidistantní). 2. Nalezení evoluční rovnice – předpis, který stanoví postup výpočtu nového stavu X(ti+1) ≡ Xi+1 v čase ti+1 z původního stavu X(ti) ≡ Xi v čase ti. Symbolicky lze zapsat tento předpis takto: Xi+1= F(Xi), kde F je obecná vektorová funkce. Posloupnost vypočtených stavů tvoří tzv. trajektorii systému a výchozí stav je dán počátečními podmínkami – proto se nalezení trajektorie nazývá počáteční úloha. Pro hmotný bod se stavovým vektorem X ≡ (x, y, vx, vy) pohybující se v rovině nalezneme evoluční rovnici snadno:

(1)

Uvedený předpis reprezentuje aproximaci pohybu mezi dvěma stavy rovnoměrným přímočarým pohybem na základě znalosti zrychlení, které můžeme odvodit z uvedeného zákonu síly. Poznámka Jedná se o tzv. Eulerovu metodu řešení počáteční úlohy, která tkví v nahrazení derivace použitím konečných diferencí. Z první složky evoluční rovnice např. plyne:

Pro hmotný bod pohybující se v homogenním gravitačním poli je podle Galileova zákona vektor zrychlení g ≡ (0, -g), máme tedy konečný tvar evoluční rovnice: Ci

16


3.1

Implementace výpočtu trajektorie

Program, který použije uvedenou evoluční rovnici pro zobrazení pohybu hmotného bodu je triviální

Program nejprve načte z listu pojmenovaného "vrh" počáteční složky stavového vektoru a velikost kroku metody:

V jednoduchém cyklu počítá z původního stavu složky nového stavu, které vzápětí ukládá do patřičných buněk v témže listu. Jako čítač řádků používáme proměnnou i. Cyklus se zastaví, jakmile obsah souřadnice y poklesne pod nulovou hodnotu (dopad na zem). Program spouštíme osvědčeným tlačítkem v listu. Po spuštění programu se automaticky doplní počátek tabulky v listu o další řádky:

17


Interaktivním způsobem můžeme vložit do listu také bodový graf, kde zobrazíme dráhu bodu jako závislost souřadnice y na souřadnici x:

3.2


Program nyní upravíme do komfortnější podoby. Nejprve změníme interaktivní list "vrh" např. takto:

18


Počáteční podmínky tedy budeme zadávat prostřednictvím velikosti rychlosti a elevačního úhlu. Data trajektorie přemístíme do přidaného dalšího listu s názvem "data". Dále přidáme pole pro výstup délky a max. dosažené výšky vrhu. Program přizpůsobíme těmto úpravám např. takto:

19


Komentář Pro výpočet složek počáteční rychlosti potřebujeme zadaný úhel převést na radiány. Zavedeme proměnnou maxY s počáteční hodnotou 0 a v cyklu výpočtu bodů trajektorie testujeme, zda souřadnice y stavu nepřekročila obsah proměnné maxY. Pokud ano, vložíme do proměnné maxY tuto novou hodnotu souřadnice. Jedná se o standardní algoritmus nalezení maxima v datové řadě. Vypočtené složky posloupnosti stavů nyní ukládáme do listu "data". Po skončení cyklu provedeme výstup poslední vypočtené hodnoty souřadnice x (délka vrhu) a obsahu maxY (max. dosažená výška). Na závěr předáme vypočtené oblasti listu "data" grafu v listu "vrh". Interaktivní list aplikace nyní vypadá takto:

3.3

Interaktivní prvek posuvník

Uvažujme experimenty, kdy bychom interaktivně hledali elevační úhel pro dosažení konkrétní vzdálenosti. V takovém případě bude výhodným prvkem formulářový posuvník (scrollbar). Do listu ho vložíme podle obrázku:

20


Posuvník je interaktivní prvek, jehož vlastnosti lze konfigurovat pomocí pravého tlačítka myši: svážeme ho s buňkou obsahující hodnotu elevačního úhlu a nastavíme mezní hodnoty v intervalu <0, 90>. Poloha táhla posuvníku má nyní okamžitý vliv na obsah buňky B3 (platí to i opačně – změna obsahu buňky vyvolá změnu polohy táhla). Prvek posuvníku navíc při změně polohy táhla generuje událost, kterou svážeme s naším programem:

Nyní se vypočítá a vykreslí nová trajektorie při každé změně v posuvníku. Můžeme tak například snadno nalézt vhodný úhel pro dosažení stanovené vzdálenosti vrhu.

21


3.4

Použití pole v počáteční úloze

Pro malý časový krok může být výpočet poněkud zdlouhavý. Program musí vyplnit někdy i tisíce řádků v listu, které následně zpracuje graf. Při změně posuvníku pak musíme počkat na odezvu i několik sekund. Tato operace by se radikálně zkrátila, kdybychom vypočtené body trajektorie ukládali do pole a předali ho do grafu tak, jako v předcházejících příkladech. Narazíme však na několik problémů: 1. Potřebnou velikost pole dopředu neznáme, počet stavů v trajektorii závisí na velikosti kroku a na počátečních podmínkách. 2. Graf zpracuje jenom velmi omezený počet bodů v poli, rozhodně nelze tímto způsobem předávat tisíce hodnot. Nezbývá než přistoupit k několika kompromisům: 1. Do pole neuložíme každý vypočtený stav – jejich hustota stejně překračuje rozlišovací schopnost grafu. V programu tedy uložíme do pole například jen každý 100. stav (pro krok 0.001s tedy zobrazíme stavy po 0.1s). 2. Před cyklem výpočtu trajektorie vytvoříme pole s dostatečně velkou dimenzí, i když některé prvky polí zůstanou nevyužity. Po skončení cyklu pole zkrátíme na skutečný počet vypočtených prvků (jinak by graf zobrazoval i nevyužité prvky). Teprve pak předáme pole k zobrazení. Navrhované úpravy znamenají, že již nemusíme používat druhý list "data". Implementace programu může být následující:

22


Komentář Obě pole musíme použit jako dynamická navzdory faktu, že stanovíme jejich dimenzi konstantou hned na začátku programu. Staticky deklarovaná pole bychom totiž nemohli později zkrátit na potřebnou délku.

23


Pro určení stého stavu použijeme standardně operátor modulo (zbytek po celočíselném dělení je nula). Do polí přidáme pro zobrazení v grafu i bod dopadu. Index polí j obsahuje po skončení cyklu maximální hodnotu, na kterou musíme zkrátit původní pole. Pro zkrácení použijeme příkaz ReDim s klíčovým slovem Preserve, kterým vynutíme zachování hodnot prvků pole při změně jeho dimenze. Jinak by se obsah pole vymazal. Nyní je aplikace tvořena jediným listem a reaguje na změnu polohy posuvníku prakticky ihned:

3.5

Balistická křivka

Šikmý vrh je samozřejmě triviální úloha, která je analyticky řešitelná (parabola). Přidáme-li k silám působícím na hmotný bod odpor prostředí, situace se podstatně změní. Analytické řešení úlohy již v takovém případě neexistuje a musíme použít numerický výpočet, avšak potřebné změny v našem programu budou minimální. Odpor prostředí se projevuje silou T(v), která působí proti směru rychlosti pohybu a je závislá na velikosti rychlosti. Podle zákona síly platí: G + T(v) = ma Obvykle předpokládáme, že síla odporu prostředí má tvar T(v) = C|v|.v, je tedy úměrná čtverci rychlosti. C je záporný koeficient odporu prostředí, který závisí na geometrii objektu a hustotě prostředí. Absolutní hodnota (velikost vektoru rychlosti) je nutná proto, aby se zachovala informace o směru rychlosti. Po dosazení za působící síly dostaneme mg + C|v|.v = ma 24


Odtud získáme velikost zrychlení pro evoluční rovnici a = g + C|v|.v / m ≡ (Cv.vx /m, Cv.vy /m – g) Samotná evoluční rovnice pro stanovení nového stavu z předcházejícího stavu bude podle vztahu (1) vypadat následovně:

Na pravé straně rovnic jsou vždy pouze složky původního stavu. Parametry úlohy se nyní rozšíří nejen o koeficient odporu prostředí, ale i o hmotnost vrženého objektu, která se v případě bez odporu prostředí vykrátí a nehraje roli (Galileův zákon). O tyto parametry musíme rozšířit vstupní list aplikace. Úprava programu spočívá pouze v přizpůsobení čtení vstupních dat a úpravu výpočtu nového stavu v cyklu (pouze tři řádky programu):

List aplikace může po úpravách vypadat např. následovně:

25


4

Simulace pohybu pružného lana

Tvar pružného lana uchyceného ve dvou bodech v homogenním gravitačním poli lze v rovnováze popsat netriviální funkcí, kterou lze za zjednodušujících podmínek analyticky vyjádřit. Zachytit dynamiku pohybu takového systému již analyticky nelze a musíme použít vhodnou numerickou metodu založenou na racionálním modelu: budeme předpokládat, že lano se skládá z nehmotných lineárních pružin připojených k hmotným bodům. Důležitým požadavkem při tvorbě modelu je jeho konvergence – při zvyšování počtu pružných lineárních úseků by se již rovnovážný tvar lana neměl od jisté kritické hodnoty dále měnit. Ukazuje se, že úlohu lze pomocí VB s využitím uživatelských možností Excelu (zejména tvorby grafu) řešit velmi efektivně.

4.1

Fyzikální model

V našem modelovém případě na každý vnitřní uzel (hmotný bod) působí gravitační síla G, pružná síla levého sousedního uzlu FL a pravého FP. Výsledná síla F pak uděluje hmotnému bodu zrychlení a. Další řešení je pak analogické počáteční úloze z předcházejícího příkladu s tím rozdílem, že pohybujících bodů je více a jejich vzájemná poloha ovlivňuje výslednou sílu působící na jednotlivé uzly – viz obrázek:

(xL,yL)

FL FP

(xP,yP)

(x,y) G

F

Úloha má zřejmě rovinný charakter, všechny vektory mají tedy pouze dvě kartézské složky. Vektory působících sil lze vyjádřit jednoduše: G = mg, kde m je hmotnost bodu, g je vektor gravitačního zrychlení (jeho orientaci zvolíme přirozeně gx = 0 a gy = -g) Velikosti pružných sil jsou: (1a) (1b) 26


kde dL a dP jsou okamžité vzdálenosti zkoumaného uzlu od levého a pravého souseda; d je rovnovážná délka pružin a k je jejich tuhost (tyto parametry předpokládáme stejné pro všechny segmenty). Okamžité vzdálenosti jsou: (2a) (2b) Pružné síly zjevně působí podél spojnice sousedních bodů, pro jejich složky podle (1a) tedy platí: (3a)

Analogicky také pro sílu pravé pružiny: (3b)

Nyní již můžeme sečíst složky všech působících sil na uvažovaný hmotný bod: (4a) (4b) Složky vektoru okamžitého zrychlení uvažovaného bodu jsou pak

,

(5)

Okamžitý stav hmotného bodu je zřejmě určen jeho polohou (x, y) a vektorem rychlosti v. Celkový vektor stavu má tedy složky (x, y, vx, vy). Pro přechod do nového stavu (x', y', v'x, v'y) po uplynutí malého časového kroku h budeme pohyb uzlu aproximovat podle jednoduché Eulerovy metody rovnoměrně zrychleným pohybem: (6)

(7)

Neustálým opakováním výpočtu nového stavu všech uzlů na základě předcházejících stavů a aktuálních působících sil simulujeme časový vývoj celého lana.

4.2

Tlumení

Takto vytvořený systém bude konzervativní – nebude konvergovat k rovnovážnému stavu. Proto obvykle přidáváme lineární tlumící sílu, jejíž vektor je úměrný okamžité rychlosti, je však opačně orientován (koeficient úměry je záporný). Vztahy (4) tedy rozšíříme na tvar: (8a) 27


(8b) Tlumicí koeficient C se tak stane dalším parametrem úlohy. Potřebnou rychlost pro výpočet výsledné síly působící na vnitřní uzel získáme ze stavového vektoru uzlu. Ostatní vztahy zůstanou zachovány.

4.3

Implementace

Jednoduchá implementace uvedené úlohy může být např. následující: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

Sub lano() ' --- deklarace proměnných ---' --- parametry úlohy Dim m As Double Dim L As Double Dim d As Double Dim h As Double Dim k As Double Dim C As Double Dim n As Integer Dim vAnim As Integer

' ' ' ' ' ' ' '

hmotnost uzlu vzdálenost bodů uchycení (délka neprověšeného lana) rovnovážná délka segmentu (nezatížené pružiny) časový krok (s) tuhost jedné pružiny koeficient tlumení počet segmentů (pružin) rychlost animace (překreslení grafu po vAnim krocích)

' --- pole stavů uzlů Dim poleX() As Double Dim poleY() As Double Dim poleVx() As Double Dim poleVy() As Double

' ' ' '

pole pole pole pole

souřadnic x uzlů souřadnic y uzlů x-ových složek vektoru rychlosti uzlů y-ových složek vektoru rychlosti uzlů

' --- čítače Dim i As Integer, j As Integer ' --- proměnné pro výpočet stavu konkrétního uzlu Dim x As Double, y As Double ' souřadnice Dim vx As Double, vy As Double ' složky rychlosti Dim ax As Double, ay As Double ' složky zrychlení Dim xL As Double, yL As Double ' souřadnice levého sousedního uzlu Dim xP As Double, yP As Double ' souřadnice pravého sousedního uzlu Dim FLx As Double, FLy As Double ' složky pružné síly levé pružiny Dim FPx As Double, FPy As Double ' složky pružné síly pravé pružiny Dim dL As Double, dP As Double ' okamžitá délka levé a pravé pružiny Const g = -9.81

' velikost gravitačního zrychlení

' --- načtení parametrů úlohy z listu --n = Worksheets("parametry").Cells(1, 2).Value m = Worksheets("parametry").Cells(2, 2).Value / (n - 1) h = Worksheets("parametry").Cells(4, 2).Value L = Worksheets("parametry").Cells(3, 2).Value d = L / n k = Worksheets("parametry").Cells(5, 2).Value / L / d C = Worksheets("parametry").Cells(6, 2).Value vAnim = Worksheets("parametry").Cells(7, 2).Value ' --ReDim ReDim ReDim ReDim

dynamická dimenze polí podle počtu segmentů --poleX(n) poleY(n) poleVx(n) poleVy(n)

28

' ' ' ' ' ' ' '

počet segmentů hmotnost uzlu časový krok délka lana rovnovážná délka pružiny tuhost segmentu tlumení rychlost animace

Řešené příklady 51 ' --- naplnění polí (počáteční podmínky) --52 For i = 0 To n 53 poleX(i) = d * i 54 poleY(i) = 0 55 poleVx(i) = 0 56 poleVy(i) = 0 57 Next i 58 59 ' --- nekonečný cyklus počítání nových stavů uzlů --60 Do While True 61 ' --- cyklus výpočtu stavů všech vnitřních uzlů --62 For i = 1 To n - 1 63 x = poleX(i) ' načtení stavu i-tého uzlu z pole 64 y = poleY(i) 65 vx = poleVx(i) 66 vy = poleVy(i) 67 68 xL = poleX(i - 1) ' souřadnice levého sousedního uzlu 69 yL = poleY(i - 1) 70 xP = poleX(i + 1) ' souřadnice pravého sousedního uzlu 71 yP = poleY(i + 1) 72 73 dL = Math.Sqr((x-xL)*(x-xL)+(y-yL)*(y-yL)) ' okamžitá délka levého segmentu 74 dP = Math.Sqr((x-xP)*(x-xP)+(y-yP)*(y-yP)) ' okamžitá délka pravého segmentu 75 FLx = k * (d - dL) * (x - xL) / dL ' vektor levé pružné sily 76 FLy = k * (d - dL) * (y - yL) / dL 77 FPx = k * (d - dP) * (x - xP) / dP ' vektor pravé pružné síly 78 FPy = k * (d - dP) * (y - yP) / dP 79 ax = (FLx + FPx - C * vx) / m ' zrychlení uzlu 80 ay = (FLy + FPy + m * g - C * vy) / m 81 82 ' --- aktualizace stavu i-tého uzlu 83 poleVx(i) = vx + ax * h ' nový vektor rychlosti 84 poleVy(i) = vy + ay * h 85 poleX(i) = x + vx * h + 0.5 * ax * h * h ' nové souřadnice 86 poleY(i) = y + vy * h + 0.5 * ay * h * h 87 Next i ' konec výpočtu jednoho časového kroku 88 89 ' --- aktualizace grafu po vAnim krocích 90 j = j + 1 91 If j > vAnim Then 92 j = 0 93 ' --- definice zdroje dat pro graf 94 Worksheets("parametry").ChartObjects(1).Chart.SeriesCollection(1).XValues = poleX 95 Worksheets("parametry").ChartObjects(1).Chart.SeriesCollection(1).Values = poleY 96 ' --- vynucení překreslení grafu 97 Application.ScreenUpdating = True 98 End If 99 100 Loop ' konec nekonečného cyklu výpočtu nových stavů uzlů 101 End Sub

29


Obr. 3.1 – Příklad listu úlohy pro zadání parametrů a výstupní animaci

4.4

Komentář

Program je velmi krátký, valnou část představují deklarace proměnných a vysvětlující komentáře. Implementace přímočaře odpovídá uvedenému fyzikálním modelu. K jednotlivým sekcím programu: Uložení stavových vektorů uzlů (řádky 14 – 18 a 45 – 49): Pro uložení stavových vektorů všech uzlů použijeme čtyři pole – pro každou složku stavového vektoru jedno. Použijeme přitom dynamická pole, neboť jejich rozsah v době deklarace neznáme. Po načtení počtu uzlů pole předimenzujeme na požadovaný počet prvků. Načtení parametrů úlohy z listu MS Excel (35 – 43): Jednoduchý formulář pro načtení parametrů vypadá např. podle obr. 3.1. Je-li n počet segmentů (pružin), pak hmotnost uzlu = (hmotnost lana) / (n – 1). Za zmínku stojí ještě tuhost jedné pružiny odvozená od celkové tuhosti soustavy. Má-li se lano chovat obdobně pro různý počet segmentů, musí být tuhost jedné pružiny odvozena od tuhosti jednotky délky tak, jak je uvedeno na řádku 41. Počáteční podmínky (51 – 57) Na počátku lano uvedeme do horizontální polohy, všechny pružiny jsou v rovnovážném stavu a uzly mají nulovou rychlost. Lano se pak vlastní vahou prověsí a pružiny protáhnou. Samozřejmě si lze představit i jinou vhodnou počáteční podmínku. Cyklus počítání nových stavů (60): Je použit nekonečný cyklus, který lze ukončit pouze kombinací kláves Ctrl– Break. Toto řešení není ideální, později si naznačíme elegantnější metodu.

30


Cyklus výpočtu nového stavu všech uzlů (61 – 87): Všimneme si, že uzly s indexem 0 a n do cyklu nezahrnujeme. Jsou to uzly, v nichž je lano upevněno, jejich souřadnice a rychlosti se nemění. Načtení proměnných pro výpočet sil v jednom uzlu (63 – 71): Podle modelu potřebujeme pro výpočet nového stavu i-tého uzlu jeho aktuální stav a souřadnice sousedů. Tyto hodnoty jsou uloženy v polích a ve výpočtu by bylo možné použít přímo prvky pole (nezavádět další proměnné x, y, xL, xP atd.). Uvedené proměnné pro jeden prvek jsou zde použity pro lepší přehlednost programu. Vypočet sil a zrychlení uzlu (73 – 80): Výpočet přesně kopíruje vztahy (1) – (8) uvedeného modelu a nevyžaduje další komentář. Výpočet nového stavu uzlu (82 – 86): Použili jsme velmi jednoduchou Eulerovu numerickou metodu, která není příliš přesná a vyžaduje tedy malý časový krok. Malá přesnost v našem případě není příliš na závadu, zajímáme se spíše o rovnovážný výsledek. Měli bychom však mít na paměti, že použitím složitější metody (např. klasická metoda Runge– Kutta) můžeme efektivitu programu snadno řádově zvýšit. Animace (89 – 98): Chceme-li sledovat vývoj našeho systému, musíme periodicky překreslovat bodový X–Y graf umístěný v témže listu, jako formulář pro vstupní parametry úlohy (viz obr. 3.1). Graf však nemůžeme aktualizovat příliš často, neboť výpočet by byl velmi pomalý. Grafický rastrový výstup je ve všech simulačních programech úzkým místem, kde procesor vyčerpá většinu přiděleného času. Proto jsme zavedli jednoduchý čítač vypočtených časových kroků a teprve po dosažení stanoveného počtu (ve vstupním listu bylo zadáno1000 kroků) graf překreslíme – objektu grafu předáme patřičná pole souřadnic x, y uzlů. Hierarchie objektů grafu byla diskutována v modulu CU01_M02, kde jsme přiřazovali řadám grafu XValues a Values patřičné rozsahy buněk v listu. Řádek 94 a 95 však naznačuje, že místo rozsahů buněk lze na toto místo vložit přímo název proměnné typu pole. Tímto způsobem náš program nemusí generovat žádná data v buňkách listů a je proto nejen jednodušší, ale zejména rychlejší. Je to jedna z mála efektivních možností pracovat s polem jako celkem. I zde si však tvůrci objektu Chart vyžádali nemalou daň – pole se ve skutečnosti předá do řady jako řetězec znaků s omezenou délkou. Zde nám tento fakt sice nevadí, ale v případě, že bychom chtěli tuto metodu použít pro velký počet uzlů, nesmíme být překvapeni ořezáním počtu bodů v grafu na relativně nízkou hodnotu. Na řádku 97 nalezneme speciální přiřazení Application.ScreenUpdating = True, které vynucuje překreslení obsahu oken aplikace. Implicitně se totiž VB chová tak, že zakáže kvůli zvýšení rychlosti překreslování po dobu běhu programu. Samotná změna dat v řadách objektu grafu v průběhu programu překreslení nevyvolá a my bychom animaci nemohli sledovat. 31


Není patrně třeba dodávat, že je zapotřebí zafixovat rozsah a jednotky v jednotlivých osách grafu. V automatickém režimu by se jinak rozsah grafu přizpůsoboval vypočteným hodnotám, což je při simulaci nevhodné:

Obr. 3.2 – Pevný rozsah os v animačním grafu

4.5

Zastavení simulačního programu

Použití násilného zastavení běhu pomocí kombinace kláves Ctrl–Break není pro uživatele příliš elegantní (vstup do ladicího režimu), navíc nelze v programu pokračovat dalším výpočtem (viz úlohy). Nejnázorněji lze přerušit simulační cyklus speciálním tlačítkem tak, jak je naznačeno na obr. 1. Stisk tlačítka je však asynchronní událostí a v simulačním programu nemáme jednoduchou možnost periodicky testovat, zda je tlačítko aktivováno. Stisku tlačítka přiřadíme tedy jiný jednoduchý podprogram, který pouze nastaví hodnotu logické proměnné konec = true. Aby byla tato hodnota přístupná i simulačnímu programu, musí být deklarována jako globální mimo těla všech programů:

Obr. 1.3 Globální proměnná a ukončovací procedura Před simulační cyklus pak přidáme příkaz konec = false a změníme podmínku cyklu na řádku 60 na testování obsahu logické proměnné konec: do while not konec 32


Pokud modifikujeme program uvedeným způsobem, zjistíme, že tlačítko nereaguje a simulační program tímto způsobem zastavit nelze. Je to způsobeno tím, že událost stisku tlačítka není zpracována – při běžícím simulačním programu opět kvůli rychlosti Excel implicitně pozastaví zpracování všech jiných událostí. Naštěstí máme možnost zpracování asynchronních událostí v programu explicitně vynutit jednoduchým příkazem DoEvents. Příkaz vyvolá zpracování fronty událostí na úkor běžícího programu, neměli bychom ho tedy volat příliš často. Vhodným umístěním pro příkaz je např. sekce překreslování grafu, která se děje v rámci simulace pouze občas:

Nyní je již možné program zastavit rozumně jednoduchým mechanizmem.

33


4.6

Úkoly

1. Prozkoumejte chování programu pro různé parametry úlohy. Všimněte si, že pro vysokou tuhost je třeba zmenšit časový krok, jinak výpočet nekonverguje. Analyzujte také vliv tlumení na stabilitu výpočtu. 2. Rozšiřte program tak, aby se po zastavení simulačního cyklu, kdy je průhyb lana ustálen, vykreslil graf deformace (protažení jednotlivých pružin). Zjistěte místo maximální deformace a odhadněte tím, v kterém místě při překročení pevnosti lano praskne. 3. Zjistěte, zda v průběhu simulace překračují hodnoty deformace hodnoty v ustáleném stavu. Analyzujte, zda dynamické chování lana může přispět k jeho prasknutí. 4. Vytvořte graf závislosti maximálního prověšení lana na jeho hmotnosti, tuhosti apod. 5. Modifikujte program tak, že v určeném uzlu (např. uprostřed) lano zatížíte dodatečnou hmotností a prozkoumejte opět průběh deformace.

34


5

Zpracování NMEA dat z měřicího přístroje

Mnoho modulů a záznamových zařízení GPS ukládá (resp. předává) naměřená data ve formátu podle standardu NMEA (námořní navigační standard). Při měření přijímací zařízení generuje množství textových řádků; na každém je uložena NMEA zpráva popsaná v dokumentaci standardu. Soubor s naměřenými daty pak vypadá následovně: $GPGGA,180756.00,4912.16613,N,01626.75199,E,1,07,1.1,452.31,M,43.8,M,,*55 $GPGSA,A,3,30,29,16,31,24,21,10,,,,,,02.2,01.1,01.9*03 $GPGSV,3,1,11,30,27,136,41,12,03,130,26,29,64,066,39,06,04,266,32*78 $GPGSV,3,2,11,16,34,304,41,31,39,237,41,24,63,086,42,21,56,181,41*77 $GPGSV,3,3,11,10,21,057,34,13,03,348,,23,02,314,*45 $GPRMC,180756.00,A,4912.1661,N,01626.7520,E,00.0,000.0,300709,03,E,A*25 $GPGGA,180757.00,4912.16611,N,01626.75192,E,1,07,1.1,452.56,M,43.8,M,,*5C $GPGSA,A,3,30,29,16,31,24,21,10,,,,,,02.2,01.1,01.9*03 $GPGSV,3,1,11,30,27,136,40,12,02,131,25,29,64,066,40,06,04,266,33*75 $GPGSV,3,2,11,16,34,304,40,31,39,237,41,24,63,086,42,21,56,181,41*76 $GPGSV,3,3,11,10,21,057,32,13,03,347,,23,02,313,*4B $GPRMC,180757.00,A,4912.1661,N,01626.7519,E,00.0,000.0,300709,03,E,A*2E $GPGGA,180758.00,4912.16610,N,01626.75188,E,1,07,1.1,452.72,M,43.8,M,,*5F $GPGSA,A,3,30,29,16,31,24,21,10,,,,,,02.2,01.1,01.9*03 $GPGSV,3,1,11,30,27,136,40,12,02,131,25,29,64,066,39,06,04,266,32*7A $GPGSV,3,2,11,16,34,304,39,23,02,313,28,31,39,237,39,24,63,086,40*70 $GPGSV,3,3,11,21,56,181,40,10,21,057,32,13,03,347,*45 $GPRMC,180758.00,A,4912.1661,N,01626.7519,E,00.0,000.0,300709,03,E,A*21 …

Obr. 5.1 – Soubor s daty NMEA Každý řádek začíná příznakem typu zprávy "$GPXXX" – některé typy obsahují údaje o horizontální poloze, jiné o vertikální poloze, některé zprávy popisují konfiguraci satelitů, z jejichž signálů byla poloha stanovena apod. Ze souboru obvykle potřebujeme extrahovat data určující čas a okamžitou polohu měřicího zařízení, z nichž pak můžeme stanovit tvar trajektorie, průběh rychlosti a další údaje o pohybu přístroje. Pro naše účely je nejvhodnější zpráva RMC (příznak $GPRMC), která obsahuje textové položky (oddělené čárkou) s následujícím významem: $GPRMC,180757.00,A,4912.1661,N,01626.7519,E,00.0,000.0,300709,03,E,A*2E Kde: RMC Recommended Min. Sentence C (minimální doporučené údaje) 180757.00 Satelitní čas měření, 18:07:57 UTC A Status, A = platný záznam 4912.1661,N 49 12.1661' – Šířka (první dva znaky stupně, ost. minuty) N Zeměpisná šířka je na severní polokouli 01626.7519,E 16 26.7519' – Délka (první tři znaky stupně, pak minuty) E Zeměpisná délka je na východní polokouli 00.0 Horiz. složka rychlosti v uzlech (1.944 m/s = 1 uzel) 000.0 Okamžitý azimut ve stupních 300811 Datum, 30. 8. 2011 03,E Úhel magnetické deklinace *6A Kontrolní suma zprávy

Obr. 5.2 – Kódování položek ve zprávě RMC Nejprve vytvoříme program, který umožní vybrat potřebný soubor s NMEA daty, vybere správné řádky a vloží do listu položky s časem měření a polohou: 35


Obr. 5.3 – Implementace programu pro načtení RMC zpráv do listu Použitý postup při čtení souboru je vysvětlen v modulu CU01_M02 v kapitole 4. Každý přečtený řádek nejprve testujeme podle úvodních znaků, zda obsahuje zprávu RMC. Pokud ano, rozdělíme řádek na položky pole poleRadek pomocí funkce Split; oddělovacím znakem je čárka (řetězcové funkce jsou diskutovány v CU01_M01, kap. 3.9). Pokud obsahuje druhá položka znak "A", je záznam validní a můžeme do patřičných sloupců a řádků v listu uložit položky pole, které vrátila funkce Split: index 1 – čas měření, 3 – zeměpisná šířka, 5 – zeměpisná délka. Položky s indexem 4 a 6 netestujeme, předpokládáme, že měření probíhá na severní a východní polokouli. Nyní si program po startu vyžádá pomocí dialogového panelu textový soubor, extrahuje z něj potřebná řádky a naplní tři sloupce listu textovými daty:

36


Obr. 5.4 – Extrahované položky z NMEA souboru Zakódovaná data musíme před dalším použitím nejprve upravit:

5.1

Satelitní čas

Řetězec času je v RMC zprávě podle obr. 5.2 zakódován šesti znaky hhmmss (resp. při vyšší frekvenci měření může následovat ještě desetinná tečka a zlomky sekund). Obsah buňky bychom měli uložit ve vnitřním formátu času, který používá VB. K tomu lze použít vestavěnou funkci TimeValue("13:45:55"), která převádí uvedené textové určení času do datového typu Date (obsahuje datum a čas). Nejprve tedy musíme vložit do textu času dvojtečky a výsledek převést zmíněnou funkcí: cas = TimeValue(Mid(RMCcas,1,2) & ":" _ & Mid(RMCcas,3,2) & ":" & Mid(RMCcas,5,2)) (použití funkce Mid – viz CU01_M01, kap. 3.9). Výhodou uložení času v datovém typu Date je zejména možnost přirozených aritmetických operací s obsahem proměnné (porovnávání časových okamžiků, přičítání, odečítání atd.).

5.2

Transformace souřadnic

Při práci s polohou obvykle pracujeme v kartézském souřadném systému, kde můžeme pohodlně zobrazit trajektorii, vypočítat rychlost, zrychlení a jiné parametry pohybu. Data GPS jsou však primárně uložena v obecném souřadném systému WGS 84 (zeměpisná délka, šířka, nadmořská výška), a proto je musíme nejprve do kartézského systému transformovat. Obecně vzato má zmíněná transformace vždy aproximativní charakter – povrch Země totiž není ani kulový natožpak rovinný. Na našem území se nejčastěji používá regionální kartézský souřadný systém S–JTSK, opírající se o síť referenčních bodů a transformace mezi WGS 84 a S–JTSK je relativně složitá. My však budeme předpokládat, že měření je omezeno na relativně malou oblast a použijeme jednoduchou lineární aproximaci jak pro vzdálenosti ve východo–západním směru a rozdíl úhlů zeměpisné délky (Δx = C1Δϕ), tak pro 37


vzdálenosti v severo–jižním směru a rozdíl úhlů zeměpisné šířky (Δy = C2Δλ). Převodní koeficienty C1 a C2 jsou vzhledem k nepravidelnosti zemského povrchu závislé na konkrétní poloze; my pro jejich určení použijeme jednoduchou aproximativní tabulku:

Obr. 5.5 – List převodních koeficientů mezi zeměpisnými úhly a vzdálenostmi Vidíme, že tabulka aproximuje zemský povrch rotačním tělesem, tj. vzdálenost dvou bodů v metrech závisí na rozdílu zeměpisných úhlů sice nelineárně, ale nezávisle na zeměpisné délce. Relativní kartézské souřadnice bodu vypočteme následujícím postupem: Údaje o souřadnicích v RMC zprávě nejprve převedeme na stupně. Podle obr. 5.2 první dva znaky vyjadřují informaci o stupních zeměpisné šířky stupne = CDbl(Mid(RMCsirka,1,2)) minuty = CDbl(Mid(RMCsirka,3,Len(RMCsirka) – 2) Ostatní znaky (jejich počet je Len(RMCsirka)-2) vyjadřují minuty a jejich zlomky. Pro extrakci patřičné části znaků jsme použili opět funkci Mid. Celkový úhel tedy bude uhel = stupne + minuty/60 U zeměpisné délky je postup analogický, pro stupně jsou v tomto případě pochopitelně použity tři znaky. Pro transformaci úhlů na rovinné souřadnice potřebujeme referenční bod (již bylo řečeno, že kartézské souřadnice lze použít pouze v omezené oblasti), který ztotožníme s počátkem souřadnic (0,0). Zvolíme oním referenčním bodem

38


první stanovený bod v souboru a pro ostatní body budeme určovat souřadnice pomocí rozdílu polohových úhlů Δϕ, Δλ určovaného a referenčního bodu.

převodní koeficient

Převodní koeficient vypočteme z tabulky pomocí lineární aproximace podle následujícího obrázku:

Ci+1 C Ci

λi

λ

λi+1

zeměpisná šířka

Obr. 5.6 – Převodní koeficientu mezi zeměpisným úhlem a vzdáleností V tabulce 5.5 nalezneme řádky, mezi nimiž leží změřená zeměpisná šířka (např. v proměnné uhel). Protože hodnoty úhlů začínají podle obr. 5.5 až od druhého řádku, bude patřičné číslo řádku dáno výrazem cisloRadku = Int(uhel)+2 (připomínáme, že funkce Int vrací celočíselnou část reálného čísla). Na tomto a následujícím řádku pak nalezneme potřebné hraniční koeficienty Ci a Ci+1. Vztah pro hledaný koeficient je pak: (rozdíl λi+1 – λi je podle tabulky 5.5 jeden stupeň). Nyní již máme všechny údaje pro výpočet kartézských souřadnic změřených poloh, takže implementaci programu z obr. 5.3 můžeme doplnit o další zpracování před ukončovací příkaz End Sub:

39


5.3

Implementace zpracování dat

Obr. 5.7 – Doplnění implementace 5.3 o zpracování naměřených dat

5.4

Komentář

Postupně doplňujeme list s přístrojovými NMEA daty o další sloupce. Nejprve v cyklu převedeme zakódovaná data na reálná čísla podle výše uvedeného postupu. Při čtení souboru jsme zjistili číslo posledního uloženého řádku v datovém listu, které je uloženo v proměnné posledniRadek a složí jako horní mez počítaného cyklu. V dalším cyklu naplníme sloupec obsahující rozdíl časů mezi dvěma po sobě následujícími měřeními. Použili jsme vestavěnou funkci DateDiff, která vrací rozdíl mezi dvěma časovými okamžiky zadanými jako parametry typu Date, jednotka rozdílu se zadává textovým parametrem ("s" – rozdíl v sekun40


dách). V cyklu vynecháváme první řádek, který pochopitelně není s čím srovnat. Přístroj měří zjevně po jedné sekundě, některá měření však jsou chybná (RMC zpráva není validní) a časový interval se pak prodlouží. V posledním cyklu počítáme z tabulky v listu "WGS84" převodní koeficienty podle obr. 5.6 pro rozdíl úhlů mezi aktuálním a referenčním (prvním) bodem. Výsledkem jsou jejich kartézské horizontální vzdálenosti v metrech. Pro výpočet absolutní vzdálenosti za sebou změřených bodů si v cyklu dočasně ukládáme hodnoty relativních souřadnic v proměnných xOld a yOld, abychom je nemuseli znovu číst z listu. Jedná se standardní algoritmický postup, který urychlí výpočet. Rychlost je vypočtena standardním postupem ze vzdálenosti po sobě následujících změřených bodů a časovým intervalem měření v = d/Δt.

Obr. 5.8 – Zpracovaná data v listu spolu s grafy rovinné trajektorie a časového průběhu rychlosti. Data byla pořízena v jedoucím vozidle na Masarykově okruhu v Brně. Trajektorie je podkreslena ortofotomapou – je zřejmé, že se s mapovým podkladem výborně shoduje. Z průběhu rychlosti vidíme značnou dynamiku pohybu, vozidlo razantně akceleruje i brzdí, na přesnost satelitního měření to však nemá vliv.

41

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ FAKULTA STAVEBNÍ INFORMATIKA II MODUL 03 ALGORITMIZACE A PROGRAMOVÁNÍ V INŽENÝRSKÝCH ÚLOHÁCH

Recommend Documents