Nejprve byla anténa QFH, pak softwarové rádio s čipsetem Realtek RTL2832U, pak byl software na ovládání rádia zvaný gqrx, který je založený na gnuradio a pak tu byl uživatel, který musel čekat, až přiletí satelit, čekat čtvrt hodiny a sehnat několik příkazů na vytvoření obrázku pomocí programu wxtoimg. Uživatele to přestalo bavit a tak si napsal skript, který sám sestaví obrázek a nahraje jej na web server. Pak ale musel napsat skript, který ty obrázky seřadí podle data pořízení. Jenže ani to mu moc nepomohlo a nezbylo nic jiného, než přemýšlet jak automatizovat úplně všechno. První problém bylo sehnat čas příletu družice na viditelnou část oblohy. Z GPredict se nějak viditelně vydolovat nedal. Díky dobré duši, která implementovala GPredict do jazyka PHP a zveřejnila kód se dala první část odškrtnout. Zdrojové kódy stáhnete zde. Predict.php všechno zařizuje, pro vaše potřeby využijte souboru visible_passes.php ve složce examples. Tento soubor si jen poupravte dle svých potřeb. Já potřeboval čas příletu a dobu, po kterou satelit zůstává na viditelné části oblohy. Tento skript neběží na webserveru, je volán cronem vždy po 4. hodině ranní, kdy začíná první část přeletů viditelných od nás a zkontroluje přelety na dalších 24 hodin. Zatím běží na klasickém PC, běh na Raspberry Pi ještě vyzkouším. Pokud skript neběží, možná jsem udělal chybu zde při úpravách, ale je téměř samovysvětlující, takže není problém jej opravit/předělat podle svých potřeb.
noaa.php
<?php
/**
* This is an example of how to use Predict for determining upcoming visible
* passes of the International Space Station. Please read the inline comments
* for details.
*
* Run this from the root checkout, not from examples, or the include paths
* will not work.
*/
require_once 'Predict.php';
require_once 'Predict/Sat.php';
require_once 'Predict/QTH.php';
require_once 'Predict/Time.php';
require_once 'Predict/TLE.php';
function sortByAos($a, $b) {
if ($a[2] == $b[2]) {
return 0;
}
return ($a[2] < $b[2]) ? -1 : 1;
}
// Track execution time of this script
$start = microtime(true);
// The observer or groundstation is called QTH in ham radio terms
$predict = new Predict();
$qth = new Predict_QTH();
$qth->alt = 300; // Altitude in meters
// Munich, example east of the meridian
$qth->lat = 48.1505; // Lat North
$qth->lon = 11.5809; // Lon East
// The iss.tle file is the first 3 lines of
// http://celestrak.com/NORAD/elements/stations.txt
// Make sure you update this content, it goes out of date within a day or two
exec("wget http://celestrak.com/NORAD/elements/weather.txt -O ~USERNAME/.wxtoimg/weather.txt");
$tleFile = file('http://celestrak.com/NORAD/elements/noaa.txt'); // Load up the NOAA from celestrak
$indexes = array(19, 18, 15);
$list = array();
foreach ($indexes as $index){
$index *= 3; //3 rows for satellite
$tle = new Predict_TLE(substr($tleFile[$index], 5, 2), $tleFile[$index+1], $tleFile[$index+2]); // Instantiate it
$sat = new Predict_Sat($tle); // Load up the satellite data
$now = Predict_Time::get_current_daynum(); // get the current time as Julian Date (daynum)
// You can modify some preferences in Predict(), the defaults are below
//
$predict->minEle = 0; // Minimum elevation for a pass
$predict->timeRes = 5; // Pass details: time resolution in seconds
$predict->numEntries = 10; // Pass details: number of entries per pass
$predict->threshold = 90; // Twilight threshold (sun must be at this lat or lower)
// Get the passes and filter visible only, takes about 4 seconds for 10 days
$results = $predict->get_passes($sat, $qth, $now, 1);
$filtered = $predict->filterVisiblePasses($results);
$zone = 'Europe/Berlin'; // Pacific time zone
$format = 'H:i mdy'; // Time format from PHP's date() function
$epoch = 'U';
// Format the output similar to the heavens-above.com website
//echo "<tr><td>Sat\t<td>AOS\t<td>ELE\t<td>LOS\t</tr>";
foreach ($filtered as $pass) {
if (round($pass->max_el) < 17){
continue;
}
$list[] = array($pass->satname, Predict_Time::daynum2readable($pass->aos, $zone, $format), Predict_Time::daynum2readable($pass->aos, $zone, $epoch), Predict_Time::daynum2readable($pass->los, $zone, $epoch), round($pass->max_el));
}
}
//delete the pass with worse max elevation from overlaping ones
usort($list, 'sortByAos');
$i=0;
while ($i < (sizeOf($list)-1)){
if ($list[$i][3] >= $list[$i+1][2]) {
if ($list[$i][4] > $list[$i+1][4]){
array_splice($list, $i+1, 1);
} else {
array_splice($list, $i, 1);
}
continue;
}
$i++;
}
//echo $wake_date > /sys/class/rtc/rtc0/wakealarm
foreach ($list as $pass) {
$cmd = "echo \"cd /home/USERNAME && ./apt ".$pass[0]. " ";
$cmd .= $pass[3] - $pass[2] - 60;
$cmd .= "\" | at -m ";
$cmd .= "\"". $pass[1] . " +1 minute\"\n";
exec($cmd);
//echo $cmd;
}
//echo "Execution time: " . number_format((microtime(true) - $start) * 1000, 2) . "ms\n";
exit;
Zde vidíte, že php skript stáhnete TLE pro sebe a wxtoimg, Predict vypočítá potřebné časy, které php zapíše do příkazu at. At spustí mnou připravený skript, který zaznamená vysílání ze satelitu a přechroupe záznam do podoby obrázku. Ten nahraje na webovou prezentaci meteostanice pomocí nástroje curl. Nezapomeňte upravit username.
Bash skript původně zpracovával nahrávky z gqrx, z této doby pochází ta skvělá konstrukce na čtení datového razítka na začátku skriptu. Využívá programu rtl_fm z balíčku od osmocom. Dále tam najdete sox (třeba nainstalovat) pro práci se zvukovým záznamem a samozřejmě WXtotImg (také doinstalovat).
apt.sh
#!/bin/bash
bandwidth=150k
login='LOGIN:PASSWORD'
#Freq +- 0.006 ?
case $1 in
15)freq=137.54M;;
17)freq=137.5M;;
18)freq=137.9065M;;
19)freq=137.092M;;
*)exit 1;;
esac
filename=`date -u +%Y%m%d-%H%M%S`.s16
timeout $2 /usr/local/bin/rtl_fm -f $freq -s $bandwidth > $filename
#timeout $2 /usr/local/bin/rtl_sdr $filename.bin -f $freq -s 1.8e6
if [ $? -ne 124 ]; then
rm $filename
exit 1
fi
if [ ! -f $filename ]; then
echo "Soubor $filename neexistuje."
exit 2
fi
#set variables
date=`echo $filename | cut -d- -f1`
time=`echo $filename | cut -d- -f2 | cut -d. -f1`
year=`echo $date | cut -c1-4`
month=`echo $date | cut -c5-6`
day=`echo $date | cut -c7-8`
hour=`echo $time | cut -c1-2`
minute=`echo $time | cut -c3-4`
second=`echo $time | cut -c5-6`
output=$year$month$day$hour$minute$second.wav
hvc=noaa$1-$year$month$day-$hour$minute$second-hvc.png
no=noaa$1-$year$month$day-$hour$minute$second-no.png
mcir=noaa$1-$year$month$day-$hour$minute$second-mcir.png
tmp=noaatmp.wav
map=map.png
thumbnail=noaa.png
#QTH
lat=48.1505
lon=11.5809
ele=300
#dawn=7
#dusk=17
#day/night is determined by python script
sox -r $bandwidth -c 1 -L $filename $tmp
sox $tmp $output rate 11025
touch -r $filename $output
/usr/local/bin/wxmap -T "NOAA $1" -L"$lat/$lon/$ele" -G /home/USERNAME/.wxtoimg -p 0 -l 0 -od "$day $month $year $hour:$minute:$second" $map
/usr/local/bin/wxtoimg -c -e MCIR -m $map -t "NOAA $1" -o $output $mcir
/usr/bin/curl -T $mcir ftp://ftp.server/ --user $login
if [ $? -eq 0 ]; then
convert $mcir -resize 25% $thumbnail
/usr/bin/curl -T $thumbnail ftp://ftp.server/ --user $login
/usr/local/bin/wxtoimg -c -e NO -m $map -A -t "NOAA $1" -o $output $no
/usr/bin/curl -T $no ftp://ftp.server/ --user $login
#hourt=`echo $hour | sed 's/^0//'`
python ./sun.py
if [ $? -eq 1 ]; then
/usr/local/bin/wxtoimg -c -e HVC -m $map -A -K -t "NOAA $1" -o $output $hvc
/usr/bin/curl -T $hvc ftp://ftp.server/ --user $login
rm $hvc
fi
rm $no $mcir $thumbnail $output
fi
rm $tmp $filename $map
exit 0
Pro určení zda je den či noc používám ne zcela přesný výpočet místního hvězdného času a jeho porovnání vypočteného času východu a západu Slunce, opět hvězdného. To vše kvůli senzoru 4, který se odesílá až hodinu po východu Slunce, zde je ještě prostor pro vylepšení.
sun.py
#!/usr/bin/python
from datetime import datetime
import time
import math
def days_between(d1, d2):
return (d2 - d1).days + (d2 - d1).seconds/84600.0
def echo_time(t):
print (math.modf(t)[1])+math.ceil((math.modf(t)[0])*60)/100
#QTH
lat = math.radians(48.1505)
lon = math.radians(11.5809)
ele = 300
utcoffset = 1
#Twilight
now = datetime.utcnow()
equinox = datetime.strptime("03-20", "%m-%d")
equinox = equinox.replace(now.year)
if days_between(equinox, now) < 0:
equinox = equinox.replace(now.year - 1)
d = days_between(equinox, now)
lamb = math.radians(360.0/365.2422 * d)
print "Lambda=", math.degrees(lamb), "d=", d
beta = 0.0
epsilon = math.radians(23.5)
delta = math.asin(math.sin(lamb)*math.sin(epsilon))
alfa = math.acos(math.cos(lamb)/math.sqrt(1-math.pow(math.sin(lamb)*math.sin(epsilon),2)))
if math.sin(lamb)*math.cos(epsilon)/math.sqrt(1-math.pow(math.sin(lamb)*math.sin(epsilon),2)) < 0:
alfa = 2*math.pi - alfa
ts = math.acos(-math.tan(lat)*math.tan(delta))
if ts>math.pi:
t1 = ts
t2 = 2.0*math.pi - ts
else:
t1 = 2.0*math.pi - ts
t2 = ts
t1 = math.degrees(t1)
t2 = math.degrees(t2)
while t1 > 360:
t1 -= 360
while t2 > 360:
t2 -= 360
print alfa
riselst = t1+math.degrees(alfa)
dawnlst = t2+math.degrees(alfa)
while riselst > 360:
riselst -= 360
while dawnlst > 360:
dawnlst -= 360
print "rise in angle",
echo_time(t1/15.0)
print "in LST",
echo_time(riselst/15.0)
print "set in angle",
echo_time(t2/15.0)
print "in LST",
echo_time(dawnlst/15.0)
#JD
r=now.year
m=now.month
d=now.day + (now.hour+(now.minute+now.second/60)/60)/24
if m > 2:
r = r - 1
m = m + 12
JD = 1720994.5 + 365.25*r + r/400.0 +- r/100.0 + 30.6*(m+1) + d + 2.0
print "JD=",JD
#JD best
unix = time.time()
J = 2440587.5 + unix/86400.0
T = (J - 2451545.0) / 36525.0
TT = 64.184 + 59.0 * T - 51.2 * math.pow(T,2) - 67.1 * math.pow(T,3) - 16.4 * math.pow(T,4)
JTT = J + TT / 86400.0
d = JTT - 2451545.0
print "JD=",JTT
print d
# Sidereal time
J2k = datetime.strptime("2000-01-01 12:00 UTC", "%Y-%m-%d %H:%M %Z")
#d = days_between(J2k, now)
print d
GMST = 280.46061837 + 360.98564736629 * d
while GMST > 360:
GMST -= 360
LMST = 280.46061837 + 360.98564736629 * d + math.degrees(lon)
while LMST > 360:
LMST -= 360
print "GMST",
echo_time(GMST/15.0)
print "LMST",
echo_time(LMST/15.0)
if LMST > riselst or LMST < dawnlst:
print "Day"
exit(1)
else:
print "Night"
exit(2)
t = GMST + math.degrees(lon) + 1.00273790935*(now.hour+now.minute/60 - 1)
print "t=",
echo_time(t/15.0)
h = math.degrees(math.asin(math.sin(delta)*math.sin(lat)+math.cos(delta)*math.cos(math.radians(LMST))*math.cos(lat)))
print "h=", h, "d=",
echo_time(math.degrees(delta))
