Wpływ sezonowych zmian w jonosferze na pomiary GPS

Od czterech lat biorę udział w pasywnym eksperymencie sejsmicznym “13BB star”. Będę na ten temat pisał jeszcze wiele razy, ale muszę tutaj zrobić dwa zdania wstępu. Celem eksperymentu jest rejestracja naturalnych zjawisk sejsmicznych, co pozwala na określanie głębokich struktur w Ziemi. W ramach eksperymentu na pomorzu pracowało 13 samodzielnych stacji poukrywanych głęboko w lasach. Każda ze stacji wyposażona była w odbiornik GPS, których celem było dostarczanie bardzo precyzyjnych informacji o czasie. Gdzieś na boku i przy okazji raz na godzinę każda stacja zapisywała pozycję według gps: długość i szerokość geograficzną oraz wysokość n.p.m.

Któregoś dnia postanowiłem wczytać do matlaba te dane dla jednej ze stacji. Gdy narysowałem wykres zależności mierzonej wysokości n.p.m. od czasu zainteresowała mnie mało widoczna zmiana sezonowa. Zupełnie nie dziwi rozrzut pomiarów o +- 10 metrów, bo pomiary “cywilnym” (jako przeciwieństwo do geodezyjnego) GPS są obarczone mniej więcej takim błędem. Zaciekawiło mnie sezonowe przesunięcie pomiarów, które jest bardzo dobrze widoczne po narysowaniu tygodniowej średniej kroczącej:

out_13bb_gps_001

Sprawdziłem od razu czy podobny efekt widać na pozostałych stacjach:

out_13bb_gps_002

Jeszcze ciekawej sprawa wygląda gdy porówna się średnie kroczące dla wszystkich stacji poprzez odjęcie średniej z całego okresu. Poniżej średnie 7-dniowe:

out_13bb_gps_003_7

30-dniowe:

out_13bb_gps_003_30

i 60-dniowe:

out_13bb_gps_003_60

Na powyższych wykresach widać wyraźną sezonowość zmian. Poprosiłem o komentarz specjalistę od pomiarów referencyjnych GPS: “jonosfera, standard, nuda”. Natomiast dla mnie ciekawe jest, że skoro efekt jest tak dobrze znany i tak powtarzalny, czemu nie została wprowadzona na to poprawka. To właśnie dynamika zmian jonosfery powoduje, że cywilne GPS mają ograniczoną dokładność. Sprzęt geodezyjny wykorzystuje stacje permanentne, dzięki którym można określi wpływ jonosfery na pomiar w danym momencie (stacja permanentna się nie rusza, więc odchylenie pomiaru wynika właśnie ze zmian w jonosferze). Takich szybkich lokalnych zmian nie da się przewidzieć, ale zmiany między latem a zimą wydają się całkowicie powtarzalne i oczywiste.

Gdyby ktoś był zainteresowany kodem użytym do wygenerowania 3 dolnych rysunków:

%% 7 dni
clf; hold on
files = dir('*.gps');

colors = lines(length(files));

for fi = 1:length(files)
	X = load(files(fi).name);
	D = X(:,1) - circshift(X(:,1),[1 0]); 
	idx = [1; find(D>1); length(X(:,1))]; % Szukanie przerw w danych większych niż 1 dzień
	X(X(:,4)>500,4)=NaN; 
	X(X(:,4)<50,4)=NaN; % Usuwanie pomiarów "od czapy"
	for a = 1:length(idx)-1
		i = idx(a):idx(a+1)-1;
		i(isnan(X(i,4)))=[];
		f = nanmean(X(:,4)); % Średnia dla całego okresu
		if length(i) > 7*24 % Badamy tylko fragmenty o długości > 7 dni
			% Obliczenie i narysowanie średniej kroczącej przesuniętej o średnią globalną
			plot(X(i,1), tsmovavg(X(i,4),'s',7*24,1)-f, 'Color', colors(fi,:), 'LineWidth',3)
		end
	end
end


% Rysowanie i podział rysunku na pory roku

axis([datenum(2013,09,01) datenum(2016,07,31), -6 4])
pory = [datenum(2013,09,22) datenum(2013,12,21) datenum(2014,03,20) datenum(2014,06,21) datenum(2014,09,23)  datenum(2014,12,21) datenum(2015,03,20) datenum(2015,06,21) datenum(2015,09,23) datenum(2015,12,22)  datenum(2016,03,20) datenum(2016,06,20) datenum(2016,09,22)];

names = {'Wiosna', 'Lato', 'Jesien', 'Zima'};
for a = 1:length(pory)-1
	text(mean(pory(a:a+1)),3.5,names{rem(a+1,4)+1},'HorizontalAlignment', 'Center', 'FontSize', 15);
	plot([pory(a) pory(a)], [-10 10],'k--')
end

set(gca,'xTick', pory)

datetick('x','YYYY-mm-dd','keepticks')

hgexport(gcf, '13BB_GPS_003_7.png', hgexport('factorystyle'), 'Format', 'png', 'resolution', 300, 'Width', 40, 'Height', 10);

| Komentarze

Burze w Polsce, edycja 2014

Pod koniec 2013 roku wpadłem na pomysł, aby regularnie zapisywać sobie obrazek (mapę) przedstawiający aktualne lokalizacje wyładowań atmosferycznych. Korzystałem wtedy ze strony blide.de. Oczywiście nie miałem na myśli robienia tego ręcznie, więc uruchomiłem na serwerze skrypt, który równo co 15 minut zapisywał aktualny obrazek. Skrypt działa do dziś i cały czas zapisuje aktualne mapy. Przypomniałem sobie o nim, gdy na serwerze zaczęło brakować miejsca na dane. Pobrałem obrazki na swój komputer i zacząłem zabawę. Na pierwszy ogień poszedł pierwszy pełny rok – 2014.

Teoretycznie w ciągu 365 dni powinienem zapisać 365*24*4 = 35040 obrazków. Faktycznie zapisałem 34999, czyli całkiem nieźle. Pierwsza, oczywista myśl to animacja. Połączyłem wszystkie obrazki w filmik:

W kolejnym korku chciałem podejść do tematu bardziej kompleksowo. Jak oszacować ilość wyładowań dla każdej z map. Oczywiście ilości się nie da, ale da się oszacować obszar objęty burzą. Na mapach kolorem czerwonym oznaczone są miejsca, gdzie wystąpiło co najmniej jedno wyładowanie w przeciągu ostatnich 15 minut. Niestety pliki są zapisane w mocno skompresowanym jpg, ale wyraźnie widać co czerwone a co nie. Napisałem na kolanie szybki skrypt w Matlabie:

clear;clc
files=dir('2014/blids/*.png');
S = zeros(length(files), 2);
for a = 1:length(files)
	name = ['2014/blids/' files(a).name];
	C = strsplit(name, '/');
	C2 = strsplit(C{end},'.');
	date = C2{1};
	time = C2{2};
	I = imread(name);
	idx = (I(:,:,1)>100 & I(:,:,2)<100 & I(:,:,3)<100);
	S(a,1) = sum(sum(idx));
	S(a,2) = datenum(str2double(date(1:4)), str2double(date(5:6)), str2double(date(7:8)), str2double(time(1:2)), str2double(time(3:4)), 0);
	
	fprintf('%d\t%d\n',a,S(a,1));
end

Dla każdej mapy skrypt obliczył i zapisał ilość pikseli, w których miały miejsca wyładowania. Teraz wystarczy to ładnie nanieść na wykres: np. sumarycznie dla każdego dnia.

clf
days = datenum(2014,01,01):datenum(2014,12,31);
storm = zeros(size(days));

for d = 1:length(days)
	storm(d) = sum(S(S(:,2)>=days(d)&S(:,2)0)=log10(storm(storm>0));
stairs(days-datenum(2014,01,01),storm, 'LineWidth',2)
axis([0 365 0 2.5e5])
hold on
for month = 1:12
	plot(ones(1,2)*datenum(2014,month,1)-datenum(2014,01,01),[0 2.5e5], 'k--')
end

A oto efekt:

001

Pionowe linie oddzielają kolejne miesiące. Widać wzmożoną aktywność burzową od kwietnia do września. Maksimum przypada na początku sierpnia, co widać na filmie.

Jak wykorzystalibyście takie dane? W wolnej chwili zrobię animację i wykres dla 2015 roku.

| Komentarze

Instalacja iOS 10 psuje iPhone!

Dziś o 19 naszego czasu Apple udostępniło nową, rewolucyjną wersję swojego systemu mobilnego oznaczoną jako iOS10. Cały świat razem ze mną rzucił się do instalacji. Zostawiłem na godzinę podłączony telefon do ładowarki, żeby pobrał i zainstalował uaktualnienie. Po powrocie zobaczyłem to:

img_0002

Dziwna sprawa, pierwszy raz zarządał podłączenia do komputera kablem, ale ok – chce to ma. Niestety po podłączeniu do komputera pojawił się komunikat, że telefon jest w trybie awaryjnym i trzeba od nowa przeprowadzić pobieranie i instalację system. Ok, jestem w stanie to przeboleć. Nie spodziewałem się jednak tego:

dsc02109

Jak żyć!?

Podobno problem dotyczył jedynie osób, które pobierały uaktualnienie zaraz po jego udostępnieniu. Żadne pocieszenie – mój telefon nie działa i nie zanosi się, żeby zadziałał w najbliższym czasie. Zawiodłem się na tobie Apple :(.

Edit (2016-09-14 03:15): Chyba się udało…

| Komentarze

Mapa Google w dużym formacie

Bardzo lubię mapy. Pewnego dnia wcześniej w tym roku wpadłem na pomysł, aby ozdobić ścianę drogową mapą Polski. W pomyśle był jeden mały haczyk: nie chciałem byle jakiej mapy, chciałem mapę Google. Szybki przegląd internetu nie przyniósł żadnych konkretnych rozwiązań – super, będę pierwszy!

Mapy Google wyświetlają różną ilość szczegółów w zależności od stopnia przybliżenia. Gdy na ekranie mieści się cała Polska, szczegółów nie ma wcale. Moja mapa ścienna miała być duża (planowałem coś w okolicach metr na metr), więc potrzebowałem szczegółów. Zdecydowałem się na napisanie aplikacji w C#, której zadaniem było wyświetlanie kolejnych fragmentów mapy i robienie screenshota. Rozwój aplikacji zakończył się w momencie, gdy zadziałała, więc nawet guziki nie były podpisane:

Dla średniej szczegółowości wystarczyło pokrycie Polski za pomocą ~120 screenów, a dla wysokiej zrobiłem ich ponad 1200.

Samo robienie screenów to łatwy fragment. Trudniej było skleić je w całość. Pamiętajmy, że ziemia jest kulą, więc samo sklejenie obrazków spowoduje powstanie dziur, ponieważ relacja pomiędzy stopniami długości geograficznej a kilometrami zależy od szerokości geograficznej. W sieci łatwo znalazłem relację pomiędzy pixelami mapy a fizycznymi odległościami w zależności od szerokości geograficznej:

Meters per pixel = 156543.03392 * Math.cos(latLng.lat() * Math.PI / 180) / Math.pow(2, zoom)

W kolejnym kroku wycinałem z moich screenów fragment mapy bez elementów interfejsu – dla ułatwienia zawsze taki sam prostokąt na środku ekranu. Robiłem to oczywiście automatycznie za pomocą polecenia convert. Przy okazji przycinania konwertowałem wycięty fragment na format tiff. Mając obliczone z wcześniejszego wzoru współrzędne geograficzne rogów fragmentu mapy konwertowałem plik graficzny tiff do formatu geotiff za pomocą polecenia gdal_translate. Geotiff to nic innego jak zwykły tiff z dodatkowymi informacjami dotyczącymi jego orientacji w przestrzeni. Ten proces był sterowany skryptem w MATLABie:

lat_min = 48.7;
lat_max = 55;
lon_min = 13.8;
lon_max = 24.6;

lat_step = .6;
lon_step = 1.2;


xof = 100;
yof = 100;
x = 1920 - 2*xof;
y = 982 - 2*yof;

for lat = lat_min:lat_step:lat_max
	for lon = lon_min:lon_step:lon_max
		im_width = 1720;
		im_height = 782;
		zoom = 10;

		meters_vertical = (156543.03392 * cos(deg2rad(lat))  / (2^zoom));

		[lat_tm, lon_tm] = reckon(lat, lon, km2deg(meters_vertical*(im_height/2)/1000), 0); 
		[lat_bm, lon_bm] = reckon(lat, lon, km2deg(meters_vertical*(im_height/2)/1000), 180); 

		meters_horizontal_top = 156543.03392 * cos(deg2rad(lat_tm)) / (2^zoom);
		meters_horizontal_bottom = 156543.03392 * cos(deg2rad(lat_bm)) / (2^zoom);

		[lat_tr, lon_tr] = reckon(lat_tm, lon_tm, km2deg(meters_horizontal_top*(im_width/2 - .5)/1000), 90); 
		[lat_tl, lon_tl] = reckon(lat_tm, lon_tm, km2deg(meters_horizontal_top*(im_width/2 - .5)/1000), 270); 

		[lat_br, lon_br] = reckon(lat_bm, lon_bm, km2deg(meters_horizontal_bottom*(im_width/2 - .5)/1000), 90); 
		[lat_bl, lon_bl] = reckon(lat_bm, lon_bm, km2deg(meters_horizontal_bottom*(im_width/2 - .5)/1000), 270); 

		system(sprintf('convert source/mapa_%8.5f_%8.5f.png -crop %dx%d+%d+%d -resize %dx%d -transparent white crop/mapa_%8.5f_%8.5f.tiff', lat, lon, x, y, xof, yof, 2*x, 2*y, lat, lon));
		system(sprintf('gdal_translate -of GTiff -a_srs WGS84 -a_ullr  %8.5f %8.5f %8.5f %8.5f crop/mapa_%8.5f_%8.5f.tiff crop/geomapa_%8.5f_%8.5f.tiff', lon_tl, lat_tl, lon_br, lat_br, lat, lon, lat, lon));
	end
end

Następnie za pomocą skryptu gdal_merge.py połączyłem wszystkie geotiffy w jednego wielkiego geotiffa. Tutaj praca mogłaby się zakończyć, ale chciałem, żeby moja mapa miała jakieś fajne odwzorowanie, ramkę itp. Najlepszym znanym mi programem do rysowania map jest pakiet GMT (http://gmt.soest.hawaii.edu/). Za pomocą prostego skryptu przygotowałem pdf (ps to prawie pdf) z mapą do wydruku:

gmt gmtset PS_MEDIA=Custom_60ix80i
gmt gmtset MAP_FRAME_WIDTH=.25i
gmt gmtset FONT_ANNOT_PRIMARY=32p,Helvetica,black
gmt gmtset COLOR_NAN=white
gmt psbasemap -R13.8/24.5/48.7/55 -Jl19.20/19.30/52/60/1:600000 -B1g.5:."": -K -P > mapa_all.ps
gmt grdimage ALL.tiff -D -R13.8/24.5/48.7/55 -Jl19.20/19.30/52/60/1:600000 -Gf -Q -P -O -K >> mapa_all.ps
gmt pscoast -R -Df -J -N1/3p,- -A0/0/1  -O -P >> mapa_all.ps

Tak wyglądał efekt końcowy po wydrukowaniu:

Nie dzielę się plikiem z mapą, bo jest bardzo duży i nie jestem pewnie czy google nie chciałoby mi urwać za to głowy.

| Komentarze (1)

Prędkość światła vs. internet

Pamiętam czasy, kiedy do internetu trzeba było dzwonić. Wszystko było piekielnie wolne i płaciło się za czas “on-line”. Moje pierwsze stałe łącze miało przepustowość 115200 kbps (kilobitów na sekundę). Dziwna liczba, ale właśnie takie było ograniczenie połączenia szeregowego między modemem a komputerem. Treści w internecie były oczywiście dostosowane do takich marnych przepustowości, ale i tak z zazdrością patrzyło się na ludzi posiadających łącze “T1”, gdzie przepustowości sięgały pojedynczych megabitów. Tamte czasy dawno minęły, coraz mniej ludzi korzysta z internetu po kablu miedzianym czy telewizyjnym, a coraz więcej ze światłowodu. Od pół roku i ja jestem szczęśliwym posiadaczem światłowodowego internetu. Tak teraz wygląda mój “speedtest”:

fiber

Sama przepustowość jest istotna gdy pobieramy dużo treści: np. strumieniowe wideo 4K czy inne dane w dużych pakietach. Na codzień inny parametr jest istotny – czas jaki potrzebuje mały pakiet danych, żeby dotrzeć od serwera do naszego komputera. Każde otwarcie strony internetowej to wymiana tysięcy takich pakietów z różnymi serwerami, ważne więc, żeby ten czas był możliwie najkrótszy. Na obrazku powyżej widać zmierzony ping na poziomie 2 ms – jest to czas od wysłania pakietu przez mój komputer do serwera, który go “odbija” do powrotu to mojego komputera. 2 ms to bardzo mało: 1/500 sekundy. Zacząłem się zastanawiać jak wygląda zależność czasu pinga (pingu?) od lokalizacji serwera. Do testu szybkości speedtest.net zawsze wybiera jeden z najbliższych serwerów, a tych na świecie są tysiące. Można oczywiście ręcznie wybrać kilka odległych, ale kto by się bawił w takie rzeczy ręcznie.

Ze strony http://www.speedtest.net/speedtest-servers.php pobrałem listę wszystkich serwerów. Wybrałem te, dla których podane były współrzędne geograficzne i rozpocząłem pingowanie. Każdego pingowałem 10 razy i wybierałem najlepszy wynik. Każdy testowany serwer przypisałem do kontynentu, naniosłem na mapę i obliczyłem odległość w linii prostej od mojego komputera do serwera. Tu pojawia się pierwszy problem – pakiet z danymi nigdy nie idzie po najkrótszej drodze, tylko po takiej, gdzie są kable. Mało tego, nigdy nie ma gwarancji, że za każdym razem pójdzie po tej samej drodze – na tym polega magia internetu.

Serwery naniosłem na mapę:

A wyniki na wykres:

Prędkość światła w próżni wynosi 299.792.458 m/s co w przybliżeniu daje 300.000 km/s. W światłowodzie prędkość ta jest mniejsza i wynosi ok. 200.000 km/s. Na wykresie narysowałem krzywe odpowiadające kilku wartościom prędkości. Widać, że prędkość przejścia pakietu do ameryki północnej wynosi około 130.000 km/s a do pozostałych kontynentów ok. 100.000 km/s. Dzieje się tak prawdopodobnie dlatego, że po drodze z europu do USA mamy ocean gdzie światłowód leży sobie prosto na dnie. W pozostały kierunkach droga jest bardziej skomplikowana i na lądzie jest dużo węzłów. No właśnie. Po drodze pakiet przechodzi przez wiele węzłów gdzie jest odpowiednio kierowany dalej – to wszystko kosztuje cenny czas!.

Nie publikuję tutaj kodu mojego programu, bo nie jest specjalnie ładny, a sam problem jest stosunkowo prosty. Gdyby jednak ktoś był zainteresowany kodem źródłowym to dajcie znać.

| Komentarze

Powrót po bardzo długiej przerwie

Moi mili czytelnicy,

IMG_0452Wstyd się przyznać, ale od ostatniego posta tutaj minęło prawie 5 lat. To gigantyczny szmat czasu, a u mnie sporo zmian. Zdążyłem się ożenić, spłodzić dwie córki i jestem o malutki kroczek od bycia doktorem (takim od nauki a nie wrzodów i innego świństwa), jak wszystko dobrze pójdzie to będę nim za kilka tygodni. Zdążyłem się również przeprowadzić itp. Ogólnie dużo zmian.

Kolejną zmianą jest moje najnowsze postanowienie o powrocie do pisania na tym blogu. Przyznam, że przez jakiś czas myślałem o nowym blogu po angielsku, ale ostatecznie postanowiłem nie porzucać tego co tutaj już zostało kiedyś stworzone. Układ jest prosty: minimum jeden post na tydzień – zawsze w środę o 20:00. Tematyka pozostaje bez zmian, będą to ciekawe projekty związanie z komputerami, nauką i techniką. Pomysłów mam bardzo dużo, sporo rzeczy jest od dawna zrobionych i tylko czekają na opisanie. W planach jest też kanał na YT, ale nie vlog – chciałbym ale nie mam czasu, dykcji i erudycji.

Zapraszam w środę o 20, będzie o prędkości światła i prędkości informacji w internecie.

| Komentarze

spamer – poziom ekspert

Na moim blogu – jak na każdym innym – można dodawać komentarze do postów. Oczywiście jest ta funkcja wykorzystywana przez spamerów do “popularyzowania” swoich stron internetowych: w zdecydowanej większości anglojęzycznych stron śmieciowych “oferujących” różnego rodzaju śmiecie (typu Rolexy po $5). Oczywiście blog wyposażony jest w system wyłapywania spamu, który nie dopuszcza do jego publikacji pod wpisami. Wpisy zakwalifikowane jako spam mogę przeglądać, żeby sprawdzić czy nie trafiły tam wartościowe komentarze (np. takie, w których ktoś przytacza wiele linków).

Treść spamerskich komentarzy bywa w różnych językach – polskim, angielskim, rosyjskim i innych. Najczęściej widać losową kolejność słów kluczowych i linków, efekty automatycznego tłumaczenia lub oklepane formułki typu “Bardzo podobał się ten wpis” (linkujące np. do wspomnianych Rolexów). To właśnie linki są dla spamerów najbardziej istotne. Na podstawie ilości linków prowadzącej do danej strony google oblicza (dokładny algorytm jest tajemnicą) rangę strony – im jest ona wyższa tym wyższa pozycja strony w wyszukiwarce.

Wystarczy wstępu, bo nie o takich spamerach chciałem dziś pisać. Wczoraj mianowicie trafiłem na komentarz zakwalifikowany jako spam, który zwrócił moją uwagę: był krótki i po polsku (poniekąd…). Brzmiał on tak:

Super post. Codziennie uczę się czegoś nowego na dziesiątkach różnych blogach. Zawsze miło czytać taki wartościowy materiał. Będę chciał wykorzystać część tego info na swoim blogu jeśli nie masz nic przeciwko. Oczywiście dodam link zwrotny do Twojego bloga. Dzięki bardzo za wszystko!

Jak się okazało komentarz został opublikowany do posta TME – rewelacja! – dziwne. Sama treść też była dla mnie dziwna – dlaczego ktoś miałby wykorzystywać to co ja przygotowuję na innym blogu? Spojrzałem więc na autora i zdziwiłem się jeszcze bardziej: przedszkole minsk mazowiecki. Potwierdziło to od razu dobrą decyzję filtra – komentarz ten był definitywnie spamem. Ponieważ komentarz przyciągnął moją uwagę przyglądałem się dalej. Został wysłany z adresu IP (94.102.57.64) należącego do firmy http://www.ecatel.net/ – holenderskiego dostawcy serwerów dedykowanych. Podany był nawet adres e-mail: przedszkolmiskmaz@gmail.pl.

Pozostała jeszcze jedna zapisana informacja o wpisie. Link do strony www – tak jak pisałem wyżej to na promocji adresów www zależy spamerom. Nie w tym przypadku, gdyż link ten prowadził do filmu na YouTube. Popychany wrodzoną ciekawością postanowiłem go sobie zobaczyć. Spodziewałem się jakiegoś śmiecia, i bardzo się zdziwiłem gdy zobaczyłem film reklamowy przedszkola:


kliknij na zdjęcie, aby przejść do filmu

Pod filmem podane są dane kontaktowe do przedszkola:

Dlaczego postanowiłem opisać ten przypadek? Nie jestem w stanie dojść, dlaczego ktoś, kto nie wie jak bardzo bezsensowne jest promowanie filmu na YouTube zadał sobie tyle trudu i spamuje za pomocą serwerów dedykowanych w innych krajach? Sądząc po ilości wyświetleń filmu na niewiele się to spamowanie zdało. Może wy macie jakieś pomysły: o co chodziło autorowi?

| Komentarze (4)

Rejestrator sejsmiczny

Tematyka rejestracji fal sejsmicznych była już poruszana przy okazji budowy i prezentacji czujnika drgań – geofonu. Jako rejestratora wykorzystałem wtedy wejście mikrofonowe karty dźwiękowej komputera. Niestety taki rejestrator ma wiele wad: ma małą czułość w stosunku do oczekiwanej amplitudy sygnału, jest zaszumiony przez otaczające go podzespoły komputera i ma małą rozdzielczość.

Profesjonalne rejestratory sejsmiczne to niezależne urządzenia pozwalające na podłączenie 1 lub więcej geofonów, które przetwarzają sygnał analogowy na cyfrowy wysokiej rozdzielczości przy zachowaniu małego poziomu szumów i perfekcyjnej synchronizacji czasu. Niestety urządzenia takie są bardzo drogie i przy cenach zaczynających się od kilku tysięcy dolarów – nieosiągalne dla typowego śmiertelnika amatora.

Od razu pojawił się pomysł (a było to co najmniej pół roku temu) samodzielnej budowy rejestratora. Projekt jest cały czas realizowany, ale najtrudniejsze już za mną. Postanowiłem wykorzystać najlepsze dostępne podzespoły – tak aby osiągnąć możliwie najwyższą jakość. Nie chcę dzielić się w tym miejscu szczegółami technicznymi – podejrzewam, że szersza publika nie jest tym bezpośrednio zainteresowana. Zawsze chętnie wysłucham Waszych pytań i w miarę swoich możliwości na nie odpowiem.

Z aspektów technicznych warto przedstawić sam układ przetwornika analogowo-cyfrowego. Wybrałem dedykowany dla sejsmiki przetwornik firmy Texas Instrument o symbolu ADS1282 – malutki układ o rozdzielczości 32 bitów (!!!). Układ ten dostępny jest tylko w obudowie SMD, która nie jest przyjazna przy prowadzeniu testów na płytce rozwojowej. Aby uporać się z tym problemem przygotowałem projekt płytki drukowanej dla tego układu i jego układów peryferyjnych, tak aby całość podłączyć do mikrokontrolera i zasilana na płytce rozwojowej. Wizualizacja płytki wyglądała następująco:


Wykonanie płytki na podstawie projektu zleciłem profesjonalnej firmie z Gdańska . Staranie i uważnie przylutowałem wszystkie elementy i gotowa do podłączenia płytka wyglądała następująco:

Gotowa płytka była dla mnie sukcesem, ale jej poprawne podłączenie do mikrokontrolera (atmega32) i odczytanie danych z przetwornika zajęło dobrych kilka dni pracy.

Starałem się modyfikować cały układ na płytce rozwojowej w taki sposób, aby na jego podstawie przygotować dużą płytkę drukowaną, gdzie z jednej strony byłoby podłączenie do geofonu a z drugiej wyjście USB do komputera. Ta płytka nie jest jeszcze gotowa, obecnie cały układ działa, ale na płytce rozwojowej. Przygotowane są nawet dwa niezależne kanały do dwóch geofonów:

Dzielę się tym z Wami ponieważ dużo pracy zostało wykonane i podstawowy cel został osiągnięty. Układ działa! Nie chcę wnikać ile pieniędzy pochłonęły te prace, bo sam bym się przestraszył. Nie udałoby się gdyby nie sprzęt pomiarowy użyczony przez Zakład Fizyki Litosfery IGF UW. Cały czas mam na uwadze dokończenie projektu (doskwiera brak czasu) i przygotowania urządzenia “na komercyjnym poziomie”, za pomocą którego możliwe będą konkretne pomiary i eksperymenty w terenie. Chciałbym aby całość była gotowa (przynajmniej kilka sztuk) do lata 2012.

Przypominam przy tym, że nie jestem elektronikiem, więc każdy z etapów prac był dla mnie pewnego rodzaju wyzwaniem i często błądzeniem we mgle…

| Komentarze (3)

Globalne ocieplenie vs. Polska

Zima trwa już dość długo więc może dla odmiany warto pomyśleć o czymś cieplejszym. Na przykład o globalnym ociepleniu. Wielokrotnie słyszeliśmy o podnoszeniu się poziomu oceanów. Niektórzy przewidują kilka metrów więcej wody niż teraz, gdy inni mówią o kilkudziesięciu czy nawet więcej.

Postanowiłem więc narysować dla Was mapę, na której widać ile terytorium Polski zostanie zalane.:


Poniżej animacja pokazująca zalewnie od 1 do 500 metrów wody:

| Komentarze

Matlab – uwaga na wydajność funkcji

Weźmy pod uwagę prosty przykład. Nasz program wymaga wielokrotnych konwersji z układu współrzędnych kartezjańskich na biegunowe.  Załóżmy, że konwersja jest dokonywana w pętli (na przykład dlatego, że każdy następny krok zależy od poprzedniego).

Matlab ma wbudowaną funkcję:

[THETA,RHO] = cart2pol(X,Y);

Do testowania wydajności przygotujemy sobie zbiór danych wejściowych:

v = randn(1000000,2)*1000;

Teraz dla każdej pary wartości (x, y) ze zmiennej v przeliczmy współrzędne:

for a = 1:length(v)
    [t r] = cart2pol(v(a,1),v(a,2));
end

Wykonanie powyższej pętli zajmuje na moim komputerze 11.692266 sekund.

Przeanalizujmy co robi funkcja cart2pol:

function [th,r,z] = cart2pol(x,y,z)
th = atan2(y,x); 
r = hypot(x,y); 
end

O ile atan2 jest operacją matematyczną, to hypot jest kolejną funkcją, która zwraca pierwiastek sumy kwadratów dwóch argumentów.

Co się zatem stanie jeśli wykonamy pętle z przykładu wyżej, ale zamiast funkcji cart2pol wywołamy obliczenia bezpośrednio?

for a = 1:length(v)
    t = atan2(v(a,2),v(a,1));
    r = sqrt(v(a,1)^2+v(a,2)^2);
end

Efekt jest identyczny. Różnica polega na tym, że czas obliczenia został zredukowany do 0.329579 sekundy, co daje przyspieszenie ponad 36 razy…

Wniosek – kod w matlabie warto analizować i optymalizować!

| Komentarze

Starsze wpisy »