0. Informacje ogólne

Polskie litery w tym pliku są kodowane zgodnie ze standardem ISO-8859-2.

Program "widoki" symuluje bieg światła w atmosferze pomiędzy dwoma obiektami.
Program "ekran" oblicza współrzędne rzutów obserwowanych obiektów
na wyimaginowanym ekranie, korzystając z wyników symulacji wykonanej przez
program "widoki".

Programy są rozpowszechniany na zasadach powszechnej licencji publicznej
GNU (GPL) w wersji 2 lub dowolnej późniejszej. Pełny tekst licencji
(oryginał engielski i nieoficjalne tłumaczenie polskie) można znaleźć pod
adresem http://gnu.org.pl.

Autor - Krzysztof Strasburger, e-mail: strasbur@chkw386.ch.pwr.wroc.pl.

Wersja 0.9.2

1. Wstęp

Obserwacje odległych obiektów, zwłaszcza gór (najlepiej z innych gór),
sprawiają niewymowną radość grupie dziwaków, do których sam mam zaszczyt
się zaliczać. Program niniejszy ma pomóc w odpowiedzi na pytanie, czy z
obiektu A, o pewnych współrzędnych geograficznych, znajdującego się na
wysokości h1 nad powierzchnią elipsoidy, opisującej w przybliżeniu kształt
Ziemi, da się zobaczyć obiekt B, który ma inne współrzędne geograficzne
i znajduje się na wysokości h2. Światło w ośrodku o zmiennym współczynniku
załamania, jakim jest atmosfera, nie rozchodzi się po liniach prostych,
dzięki czemu w sprzyjających warunkach można nawet zajrzeć poza horyzont.
Dokładniej, przyczyną odchylenia swiatła jest niezerowy gradient wspólczynnika
załamania światła (n) w kierunku prostopadłym do kierunku rozchodzenia się
fali świetlnej.

W programie wykorzystano następujace modele fizyczne:
- zasada Huyghensa (lub zasada Fermata - wnioski są matematycznie
  równoważne) do obliczania ugięcia promieni światła;
- model Lorentza polaryzacji dielektryka, do obliczania n w funkcji
  koncentracji cząsteczek gazu;
- równanie stanu gazu doskonałego, połączone z modelem Lorentza, daje
  zależność n od ciśnienia i temperatury;
- wzór barometryczny, opisujący ciśnienie atmosferyczne na danej wysokości.

Kluczowymi parametrami w obliczeniach są: pionowy rozkład temperatury
i ciśnienie na danej (jednej) wysokości.

2. Kompilacja programu

Program jest rozprowadzany w postaci kodu źródłowego. Do kompilacji
potrzebny jest kompilator języka C i dostępna we wszystkich systemach
standardowa biblioteka libc. W niektórych środowiskach może być potrzebna
biblioteka funkcji matematycznych libm. Zalecany jest program GNU make do
automatycznej kompilacji. Przy odrobinie szczęścia (czyli w systemach
uniksowych, w których jest zainstalowany kompilator gcc) wystarczy wydać
polecenie "make". Jeśli dostępny jest inny kompilator, może być konieczna
edycja pliku Makefile. Zresztą, zawsze warto do tego pliku zajrzeć i dostosować
parametry kompilacji do procesora własnego komputera i wersji kompilatora.
Jeśli programu "make" nie ma w systemie, można w ostateczności kompilować
źródła ręcznie, z linii poleceń, choć w miarę zwiększania liczby plików
źródłowych może to być coraz bardziej kłopotliwe.

Program został napisany i przetestowany pod linuksem, w innych systemach
mogą się pojawić niespodzianki. Użytkowników tych innych systemów proszę o
informacje, gdyby takowe rzeczywiście wystąpiły. a jeszcze lepiej - o
wskazanie miejsca w programie i przyczyny kłopotów.

Użytkownikom systemów rodziny Windows, nie mającym zainstalowanego w
systemie kompilatora języka C, zalecam pakiet mingw. Nie, nie wypróbowałem
go, bo nie używam żadnego z systemów rodziny Windows, ale jest to port
kompilatora gcc i innych programów, w tym niezbędnego do kompilacji
kodu źródłowego pakietu binutils i programu make. Powiodła się natomiast
kompilacja za pomocą kross-kompilatora mingw działającego pod linuksem,
otrzymany program działa w środowisku wine i dostałem od kolegi sygnał, 
że podobnie jest w systemie windows (nie wiem, których, alo powinno działać
we wszystkich).

3. Przygotowanie danych i obliczenia

3a. Program "widoki"

Dane są czytane ze standardowego wejścia, wyniki - wypisywane na standardowym
wyjściu. Najwygodniej jest przygotować plik z danymi i przekierować
standardowe deskryptory plików, używając polecenia:
widoki [argumenty programu] < dane > wyniki.

Argumenty programu (podawane w linii poleceń):
-h - krótki opis argumentów programu i zakonczenie dzialania,
-g - rezygnacja z obliczeń czysto geometrycznych (biegu światła po prostej),
-c - ciche obliczenia; nie są wypisywane wyniki posrednie
-d - wypisanie wielu komunikatów z obliczeń; używane w razie podejrzeń
     o wystąpienie błędów (debug).

Opis poleceń, których można używać w danych programu (czytanych ze
standardowego wejścia):

# oznacza początek komentarza, ciągnącego się do końca linii

e półoś_równikowa półoś_biegunowa 
   Wczytanie własnych parametrów elipsoidy. Domyślnie używana jest
   elipsoida Krasowskiego.

p ciśnienie_na_poziomie_morza
   Wartość w hektopaskalach, jak w prognozie pogody ;), domyślnie 1013.25.

t 1 t0 dt
   Pierwszy model pionowego rozkładu temperatury; t0 oznacza temperaturę na
poziomie morza (w stopniach Celsjusza), dt - zmianę temperatury przy
zwiększeniu wysokości o 1 m (pionowy gradient); domyślnie 0 i -0.0065.

t 2 dt1 n dt2
   Drugi model pionowego rozkładu temperatury; dt1 oznacza pionowy gradient
temperatury w dolnej warstwie powietrza, dt2 - w górnej warstwie; pomiędzy
nimi podajemy (za powyższą linią) n par wartości:
 wysokość    temperatura
przy czym wysokości w kolejnych punktach muszą być coraz większe.
W górnej i dolnej warstwie powietrza temperatura w tym modelu zmienia się
liniowo, a w warstwie środkowej jest opisywana wielomianami sklejanymi 3-go
stopnia. W liniach danych opisujących rozkład temperatury nie można używać
komentarzy (znaku #).

h d0 a0 hmax
   Przesunięcie profilu temperatury w pionie, w zależności od odległości
od punktu obserwacji (punkt ten musi już być zdefiniowany); zależność
T(h) jest zastępowana przez T(h+dh), gdzie dh=hmax*exp(-((d-d0)/a0)^2));
d jest odległością punktu bieżącego od punktu obserwacji. Wartości
parametrów d0 i a0 należy podać w kilometrach, hmax - w metrach.

w szerokość_geograficzna długość_geograficzna h1 h2 dh
  Wykonanie obliczeń ciśnienia, temperatury i ich pionowych gradientów w
  zakresie wysokości od h1 do h2, z krokiem dh, dla punktu o danych
  współrzędnych geograficznych.

s długość_pojedynczego_kroku_symulacji
  Wartość w metrach, domyślnie 1000. Co tyle metrów drukowane są pośrednie
  wyniki obliczeń.

m szerokość_geograficzna długość_geograficzna wysokość nazwa
  Miejsce w terenie. Pierwsze podane miejsce staje się punktem obserwacji,
  przy podawaniu następnych wykonywane są symulacje przebiegu światła
  z miejsca obserwacji do docelowego (w świecie realnym jest na odwrót,
  ale to nic nie zmienia w wynikach).

k
  Zakończenie wykonywania programu.

Przykładowe pliki wejściowe, opisujące różne mniej lub bardziej szalone
widoki, z komentarzami, znajdują się w katalogu "dane".

Domyślnie program wykonuje najpierw obliczenia dla światła rozchodzącego się
po prostej i i podaje odległości od miejsca obserwacji, wspólrzędne
geograficzne i wysokości, na których biegnie w tym miejscu swiatło. Potem
następuje symulacja ugięcia światła - przeważnie 3-4 przebiegi (bo światło
musi "trafić" w punkt docelowy), podczas których drukowane są te same
informacje. W ten sposób (dopiero w ostatnim przebiegu, ale przed jego
wykonaniem nie było wiadomo, że to już ostatni) można prześledzić bieg
światła i sprawdzić na mapie, czy znajdują się po drodze przeszkody
na tyle wysokie, by uniemożliwić obserwację.
Na końcu symulacji drukowane są: odległość pomiędzy punktami, najmniejsza
znaleziona wysokość, na której biegło światło, odległość od punktu
obserwacji, na której osiągnięto powyższe minimum i 3 składowe wektora,
wzdłuż którego światło wyszło z punktu obserwacji w ostatnim przebiegu
symulacji. Jeżeli użyty zostanie argument "-c", to zostanie podana tylko
ta końcowa informacja. Ma ona w zamyśle posłużyć do konstruowania panoram.
Odpowiedni program jest teraz opracowywany i powinien znaleźć się w
następnym wydaniu.

Program nie ma nadmiernie rozbudowanej diagnostyki błędów, zwłaszcza przy
wprowadzaniu danych. Użytkownik ponosi wszelkie konsekwencje ;).

3b. Program "ekran"

Program "ekran" wykorzystuje wyniki symulacji jako dane wejściowe i oblicza
współrzędne obserwowanych obiektów na płaszczyźnie (ekranie) tak, jak by je
widział obserwator. Pierwsze dwa punkty wyznaczają początek układu
współrzędnych na ekranie i dodatnią półoś y (linię pionu). Początek układu
wyznacza jednocześnie linię horyzontu. W celu wykonania poprawnych (tzn.
nadających się do wykorzystania w ewentualnej późniejszej wizualizacji)
obliczeń programem "ekran", należy najpierw dobrać punkty do symulacji
w programie "widoki".

Na ten temat należy napisać kilka słów komentarza.

W programie "widoki" nie ma żadnych ograniczeń na rozchodzenie się światła.
Ośrodek jest idealnie przezroczysty i nie ma w nim przeszkód terenowych. Te
ostatnie można stosunkowo łatwo opisać, wykonując taką samą symulację dla
obiektu, który nas interesuje i wszystkich tych, które mogłyby go ewentualnie
zasłonić. Obrazy tych przeszkód znajdą się albo wyżej nad linią horyzontu,
niż obraz obserwowanego obiektu - wtedy go zasłonią, albo niżej - wtedy nasz
obiekt powinien być widoczny. Wreszcie, obraz obiektu może się znaleźć pod
linią horyzontu - wtedy nie będzie widoczny, choćby żadna inna przeszkoda
go nie zasłaniała. Tylko właśnie - na jakiej wysokości jest linia horyzontu?
Nie symulujemy obserwacji nad oceanem i na ogół między obserwowanym obiektem
a obserwatorem nie rozciągają się bezkresne równiny, a w dodatku powietrze
nie jest idealnie przezroczyste. Zwłaszcza w warunkach inwersji wszelkie
mgły i pyły gromadzą się w nisko położonej warstwie powietrza. Możemy
wreszcie prowadzić obserwację ponad morzem chmur (którego występowanie też
zresztą wiąże się z inwersją). Gdyby w programie było uwzględnione
rozpraszanie światła, słaba przezroczystość tych warstw powietrza
wychodziłaby w sposób naturalny, w postaci zmniejszania natężenia światła
docierającego do obserwatora. Tyle, że tego zjawiska w programie nie ma,
przynajmniej na razie... Krótko mówiąc, musimy się umówić, że poniżej pewnej
wysokości nad poziomem morza światło w powietrzu się nie rozchodzi, tzn. jest
całkowicie rozpraszane lub pochłaniane. Powinniśmy dobrać tę wysokość do
warunków atmosferycznych. Program "widoki" podaje najmniejszą wysokość,
na której po drodze od obserwatora do punktu obserwowanego rozchodziło się
światło. Jeżeli jest mniejsza od przyjętej wysokości granicznej, należy
przypuszczać, że obiekt nie jest widoczny. Na wysokości granicznej będzie
się zatem znajdować umowna linia horyzontu.

Wyobraźmy sobie teraz, że symulujemy widzialność jakiejś grupy górskiej.
Po drodze z miejsca obserwacji znajdują się przeszkody. Wykonujemy symulacje
rozchodzenia się światła dla możliwie wielu punktów, mozolnie odrzucamy te
spośród nich, dla których światło zeszło po drodze poniżej ustalonej wysokości
granicznej, a co do pozostałych - musimy jeszcze ustalić, które punkty są
zasłonięte przez inne. Niewdzięczna robota. Z myślą o ułatwieniu tego zadania,
program "widoki" podaje wektory styczne w punkcie obserwacji do promieni
światła prowadzących do obserwowanych obiektów. Wektory te pochodzą
z ostatniego przebiegu symulacji, w której promień światła "trafia"
w obserwowany punkt. Program "ekran" przelicza te wektory na współrzędne
punktów na płaskim, wyimaginowanym ekranie, znajdującym się w odległości 1 km
(albo innej, bo może ona być zmieniona przez użytkownika) od obserwatora.
Dwa dodatkowe wektory wyznaczają początek związanego z ekranem,
prawoskrętnego, prostokątnego układu współrzędnych i dodatnią półoś y,
czyli kierunek pionowy. Pierwszy z nich jest obliczany jako średnia
wszystkich wektorów wskazujących obrazy obserwowanych obiektów, drugi jest
na razie wybierany trochę "na pałę", w związku z czym trzeba go później
korygować. Program "ekran" umożliwia zmiany położenia i orientacji ekranu,
według poleceń użytkownika. Każda taka zmiana pociąga za sobą zmiany
współrzędnych rzutów obrazów obserwowanych punktów na ekran.

Wyniki symulacji muszą zostać zapisane w pliku:
widoki -c -g < dane > wyniki_symulacji
po czym należy wczytać je do programu ekran, poleceniem:
ekran wyniki_symulacji
i można poddawać obróbce obrazy punktów na ekranie.
Program oblicza wektor prowadzący do centrum obrazu (tam będzie początek
układu współrzędnych związany z ekranem) i jakiś wektor pionu (na ogół
niewiele wart, ale to można zmienić). Następnie wyświetlany jest znak
zachęty (ekran>) i interpretowane są polecenia wpisywane przez użytkownika:
t tryb_prezentacji - wybór trybu prezentacji, być może za pomocą jakiego
  programu (backend); na razie "term" (wbudowany) lub "gnuplot"; może
  wymagać dodaktowych argumentów (patrz niżej dla trybu "gnuplot",
rc kąt - obrót ekranu o podany kąt wokół wektora prowadzącego do środka
  ekranu (odpowiada obrotowi aparatu fotograficznego),
rv kąt - obrót wokół osi pionowej (zmiana kierunku patrzenia),
rh kąt - obrót wokół osi poziomej (pochylanie aparatu),
d odległość - zmiana odległości ekranu od obserwatora,
wc x y z - wczytanie współrzednych wektora wskazującego środek ekranu
wp x y z - wczytanie wektora pokazującego na ekranie punkt, przez który
  przechodzi wychodzący ze środka ekranu pion,
p - pokazanie punktów (program w trybie "term" wyświetla na ekranie ich
  współrzędne, odległości od obserwatora, wysokości nad poziomem morza,
  najmniejsze wysokości biegu światła na drodze do obserwatora i nazwy),
s nazwa_pliku - zachowaj bieżące wyniki, podawane przez polecenie p, w pliku
h lista nazwa obiektów - wymienione obiekty nie mają być pokazywane
u lista nazwa obiektów - wymienione obiekty mają być pokazywane
l obiekt dy [dx] - przesunięcie podpisu obiektu w pionie i w poziomie
  względem jego bieżącej pozycji
c polecenie - komenda przekazywana do programu zewnętrznego
k - koniec obliczeń.

Program ekran w trybie "term" wyświetla współrzędne punktów jako kolumny liczb.
Prezentacja wizualna jest możliwa dzięki trybowi "gnuplot". Program ekran
współpracuje z programem gnuplot (należy go zatem najpierw zainstalować
w systemie - szczegóły na http://www.gnuplot.info). Najprostszy sposób pracy
wymaga podania nazwy roboczego pliku, służącego do zapisywania danych
przeznaczonych do wykreślenia, jako argumentu przy wyborze trybu:
t gnuplot plik_roboczy
co umożliwia wyświetlanie obrazów obserwowanych punktów w prostokątnym układzie
współrzędnych, za pomocą polecenia "p". Ponadto akceptowane jest polecenie:
c komenda_programu_gnuplot - przekazujące podaną komendę do programu gnuplot.

Jeżeli użytkownik ma zainstalowany w systemie program gnuplot w wersji
przynajmniej 4.2, może zaznaczyć punkty na tle zdjęcia, z którym zostaje
związany układ współrzędnych. Zdjęcie należy jednak najpierw przetłumaczyć
do formatu surowej mapy punktów - gnuplot nie jest przeglądarką plików
graficznych i nie "rozumie" formatów tych plików, przynajmniej na razie.
Np. plik w formacie jpg należy przetłumaczyć na format pnm:
djpeg -pnm plik.jpg > plik.ppm
a potem obciąć nagłówek pliku:
dd if=plik.ppm of=plik.rgb bs=1 skip=rozmiar_naglowka
przy czym rozmiar nagłowka można obliczyć nastepująco:
rozmiar_nagłówka=rozmiar_pliku-nx*ny*3
gdzie nx i ny oznaczają rozdzielczości w poziomie i pionie. Nagłówek pliku
ppm może mieć różne rozmiary, np. program djpeg daje 15 bajtów, natomiast
gimp - 60 lub 61 bajtów, stąd konieczność obliczenia przed dalszymi operacjami
na pliku. Tryb prezentacji należy wybrać z następującymi argumentami:
t gnuplot plik_roboczy nx ny.
Początek układu współrzędnych zostanie umieszczony w środku zdjęcia. Jeśli z
jakiegoś powodu ma być inaczej (np. posługujemy się wyciętym fragmentem
zdjęcia), można podać współrzędne lewego dolnego rogu (x0, y0):
t gnuplot plik_roboczy nx ny x0 y0.

4. Ograniczenia, niedociągnięcia i plany na przyszłość

Ciśnienie i temperatura zmieniają się tylko w pionie. Wprowadzenie poziomych
składowych gradientów tych wielkości mogłoby być interesujące, ale wymaga
wpierw odrobiny pracy ołówkiem na kartce papieru. Gradienty te mogą
mieć różne kierunki, a to trochę skomplikuje obliczenia. Będzie można wtedy
np. postawić pytanie: A co będzie, jeśli po tym zboczu, które zasłania
widok, spłynie masa zimnego powietrza? Na etapie koncepcji, tym bardziej, że
o odpowiednie dane do obliczeń wcale nie jest łatwo.

Dane pomiarowe dotyczące pionowego rozkładu temperatury można znaleźć w
internecie (http://www.uwyo.edu/upper_air) i już zostały
wykorzystane w symulacjach (np. plik slawkowski-borzawa.dat).
Dodana została możliwość przesuwania pionowego profilu temperatury
w zależności od odległości od punktu obserwacji (model pozostaje w zasadzie
jednowymiarowy). Zmiana ta powoduje pewne nieścisłości w obliczaniu
ciśnienia na danej wysokości.

Fizyka problemu - można pokusić się o uwzględnienie rozpraszania Rayleigha.
Rachunkowo nie jest to specjalnie trudne. Sztuką jest przetłumaczenie
informacji o zapyleniu i wilgotności powietrza na wartości parametru
rozpraszania. Możliwe byłoby oszacowywanie zasięgu widzialności w różnych
warunkach.

Składnia poleceń programu ekran powinna zostać porządnie zaprojektowana.
Na miejscu byłoby wykorzystanie biblioteki readline, tak by możliwa była
edycja i historia komend (teraz jest używana funkcja fgets biblioteki libc).
Zresztą, program ten wyrósł z potrzeby chwili i jest znacznie mniej
dojrzały od części obliczeniowej.

5. Podziękowania

Opisywany tu program został zainspirowany wcześniejszą o mniej więcej rok
pracą Andrzeja Wojciechowskiego, mojego starszego kolegi i partnera z gór.
Andrzej napisał program, w którym wykorzystywane są tablice refrakcji,
opracowane w oparciu o wyniki pomiarów dokonanych w prawdziwym świecie.
Dziękuję mu za udostępnienie wyników swoich obliczeń, dzięki czemu mogłem
zweryfikować założenia modelu fizycznego, leżącego u podstaw mojego kodu.

Dziękuję Mateuszowi Dłutko z zespołu Dalekich Obserwacji za zgodę na
udostępnienie programu na stronie WWW http://www.dalekieobserwacje.eu.

Dziękuję użytkownikowi Wicio z forum Dalekich Obserwacji za wskazanie
dostępnego w internecie źródła danych z pomiarów sondowania atmosfery.

Życzę wszystkim użytkownikom miłego symulowania i zapraszam do współpracy
przy rozwijaniu programu.
                                            Krzysztof Strasburger