[tracking] modernize satfind

* use argparse
* use locations.ini
* estimate doppler shift
* make visibility angle configurable

tracking/satfind3.py

@@ -1,104 +1,144 @@
 #!/usr/bin/env python3
 import sys
 import math
-import time
 from datetime import datetime, timezone, timedelta
-import getopt
+import argparse
+import re
 import os
 import termios
 import select
-from skyfield.api import load, utc, Topos
+from skyfield.api import load, Topos
 from numpy import identity
+from configparser import ConfigParser
+import pyproj
+
+ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
+lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
+
+to_lla = pyproj.Transformer.from_proj(ecef, lla)
+to_ecef = pyproj.Transformer.from_proj(lla, ecef)
 
 degrees_per_radian = 180.0 / math.pi
+speed_of_light = 299792458
+
+args = None
+
+def parseangle(value):
+    z = re.match(r"(?P<sign>[+-]?)(?P<deg>\d+(\.\d+)?)(°((?P<min>\d+(\.\d+)?)')?((?P<sec>\d+(\.\d+)?)\")?(?P<dir>[NESW]?))?$", value)
+#    print(value,end=" -> ")
+    if z is None:
+        raise ValueError("could not convert string to angle: '%s'"%value)
+    parsed = z.groupdict()
+    result = float(parsed['deg'])
+    if parsed['min'] is not None:
+        result += float(parsed['min'])/60
+    if parsed['sec'] is not None:
+        result += float(parsed['sec'])/60/60
+    if parsed['dir'] == 'S' or parsed['dir'] == 'W' or parsed['sign'] == '-':
+        result = -result
+    return result
+
+assert parseangle("123") == 123
+assert parseangle("-13.30") == -13.30
+assert parseangle("13°") == 13
+assert parseangle("13°30'S") == -13.5
+assert parseangle("13°30'S") == -13.5
+assert parseangle("13°30'1\"") == 13 + 30/60 + 1/3600
+
+def read_observer(location):
+    observer = {}
 
-options, remainder = getopt.getopt(sys.argv[1:], 'vhm:r', [
-                                                         'pos=',
-                                                         'tle=',
-                                                         'minutes=',
-                                                         ])
-
-pos='ata'
-tlefile='iridium-NEXT.txt'
-minutes=60
-verbose=False
-items=3
-repeat=False
-
-for opt, arg in options:
-    if opt in ['--tle']:
-        tlefile=arg
-    elif opt in ['--pos']:
-        pos=arg
-    elif opt in ['-m', '--minutes']:
-        minutes=float(arg)
-    elif opt in ['-r']:
-        repeat=True
-    elif opt in ['-v']:
-        verbose=True
+    config = ConfigParser(converters={'angle': parseangle})
+    config.read(['locations.ini', os.path.join(os.path.dirname(__file__), 'locations.ini'), os.path.join(os.path.dirname(__file__), '..', 'locations.ini')])
+
+    if location not in config:
+        print("Location %s not defined" %location, file=sys.stderr)
+        print("Available locations: ", ", ".join(config.sections()), file=sys.stderr)
+        sys.exit(1)
+
+    if 'name' in config[location]:
+        observer['name'] = config[location]['name']
+        args.loc = observer['name']
     else:
-        raise Exception("unknown argument: "+opt+".")
+        observer['name'] = location
+
+    if 'lat' in config[location]:
+        lat = config.getangle(location, 'lat')
+        lon = config.getangle(location, 'lon')
+        alt = config.getfloat(location, 'alt')
+        observer.update(lat=lat, lon=lon, alt=alt)
+
+#        x, y, z = to_ecef.transform(lon, lat, alt, radians=False)
+#        observer['xyz'] = np.array([x, y, z])/1000
+    elif 'x' in config[location]:
+        x = config.getfloat(location, 'x')
+        y = config.getfloat(location, 'y')
+        z = config.getfloat(location, 'z')
+#        observer['xyz'] = np.array([x, y, z])/1000
+
+        lon, lat, alt = to_lla.transform(x, y, z, radians=False)
+        observer.update(lat=lat, lon=lon, alt=alt)
+    else:
+        print("Location %s has no location information" %location, file=sys.stderr)
+        sys.exit(1)
+
+    return observer
+
+def get_locations():
+    config = ConfigParser()
+    config.read(['locations.ini', os.path.join(os.path.dirname(__file__), 'locations.ini'), os.path.join(os.path.dirname(__file__), '..', 'locations.ini')])
+
+    if config.sections():
+        return config.sections()
+
+    raise SystemExit("locations.ini missing or empty")
+
+def parse_args():
+    global args
+
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument("-v", "--verbose",   action="store_true",                        help="verbose output")
+    parser.add_argument(      "--debug",     action="store_true",                        help="debug output")
+    parser.add_argument("-l", "--loc",       choices=get_locations(), default="default", help="observer location")
+    parser.add_argument("-m", "--minutes",   type=int, default=60,                       help="minutes in the future to scan")
+    parser.add_argument("-r", "--repeat",    action="store_true",                        help="interactive mode")
+    parser.add_argument("-a", "--angle",     type=float, default=18,                     help="elevation angle")
+    parser.add_argument("-t", "--tlefile",   default="iridium-NEXT.txt",                 help="TLE filename")
+
+    parser.add_argument("items", default=3, nargs='?', type=int, help=argparse.SUPPRESS)
+
+    args = parser.parse_args()
+
+    return args
 
-if len(remainder)>0:
-    items=int(remainder[0])
 
 def loadTLE(filename):
     satlist = load.tle_file(filename)
-    if verbose:
+    if args.verbose:
         print("%i satellites loaded into list"%len(satlist))
     return satlist
 
 def print_sat(date, sat):
     print('%s %20s: altitude %4.1f deg, azimuth %5.1f deg, range %5.1f km' %
         (date, sat.name, sat.alt * degrees_per_radian, sat.az * degrees_per_radian, sat.range/1000.))
 
-satlist = loadTLE(tlefile)
-ts=load.timescale(builtin=True)
-
-def get_pos(pos):
-    locations={
-        'club':  { 'lat': 48.153543,            'lon': 11.560702,                'alt': 516 },
-        'ata':   { 'lat': 40+49./60+ 3  /60/60, 'lon': -(121+28./60+24  /60/60), 'alt': 1008},
-        'ata1h': { 'lat': 40+48./60+59.1/60/60, 'lon': -(121+28./60+18.6/60/60), 'alt': 1008},
-    }
-
-    if pos not in locations:
-        print("Location",pos,"unknown. Known locations: ",",".join([str(x) for x in locations]))
-        sys.exit(1)
-
-    lat=locations[pos]['lat']
-    lon=locations[pos]['lon']
-    alt=locations[pos]['alt']
-
-    return Topos(latitude_degrees=lat, longitude_degrees=lon, elevation_m=alt)
-
-recpos=get_pos(pos)
 
+timescale=load.timescale(builtin=True)
 def mytsutc(tm):
     """ speed up skyfield a bit / https://github.com/skyfielders/python-skyfield/issues/389 """
-    t=ts.utc(tm)
+    t=timescale.utc(tm)
     t.gast = t.tt * 0.0
     t.M = t.MT = identity(3)
     return t
 
-now=datetime.now(timezone.utc)
-tnow = mytsutc(now)
-
-sat0=satlist[0]
-days = tnow - sat0.epoch
-if verbose:
-    print('TLE file is %.2f days old'%days)
-
-if abs(days)>5:
-    print('SAT EPOCH too far away. (%.2f days)'%days, file=sys.stderr)
-    sys.exit(-1)
 
 def generate_events(satlist, tnow=mytsutc(datetime.now(timezone.utc))):
     t0 = mytsutc(tnow.utc_datetime() - timedelta(minutes=5))
-    t1 = mytsutc(tnow.utc_datetime() + timedelta(minutes=minutes))
+    t1 = mytsutc(tnow.utc_datetime() + timedelta(minutes=args.minutes))
     satev=[]
     for sat in satlist:
-        t, ev = sat.find_events(recpos, t0, t1, altitude_degrees=18.0)
+        t, ev = sat.find_events(recpos, t0, t1, altitude_degrees=args.angle)
         if len(t)>0:
             t2, ev2 = sat.find_events(recpos, t0, t1, altitude_degrees=0)
             events=[]
@@ -127,7 +167,6 @@ def deltatstr(ts):
         out="%5s ago"%deltat(now,ts)
     return out
 
-
 def format_events(satev, nitems=3):
     ad=30 # azimuth-maximum delta for north/south/east/west
     nice=[]
@@ -140,13 +179,13 @@ def format_events(satev, nitems=3):
         output=[]
         for ti, event in item[2]:
             add=""
-            name = ('rise', 'above 18°', 'culminate', 'below 18°', 'set')[event]
+            name = ('rise', 'above %d°'%args.angle, 'culminate', 'below %d°'%args.angle, 'set')[event]
             name="%-10s"%name
             if event==1:
                 sts=ti.utc_datetime()
             if event==3 and sts is not None:
                 add+=" [visible for %s]"%deltat(ti.utc_datetime(),sts)
-            if verbose and (event==3 or event==5):
+            if args.verbose and (event==3 or event==5):
                 (el,az,dist)= (item[1]-recpos).at(ti).altaz()
                 add+=" [az=%3d°]"%az.degrees
             if event==0:
@@ -168,24 +207,48 @@ def format_events(satev, nitems=3):
     return nice
 
 def print_events(nice):
-    global now
     t_=mytsutc(now)
     for sat,out in nice:
         print()
         print("%s:"%sat.name, end="")
-        (el,az,dist)= (sat-recpos).at(t_).altaz()
-        print(" [az=%5.1f° el=%5.1f° dist=%5dkm/%4.1fms]"%(az.degrees,el.degrees,dist.km,dist.m/ 299792.458),end="")
+#        (el,az,dist)= (sat-recpos).at(t_).altaz()
+        (el,az,dist,el_r,az_r,dist_r)=(sat-recpos).at(t_).frame_latlon_and_rates(recpos)
+        shift=(speed_of_light-dist_r.m_per_s)/speed_of_light*1.626e9-1.626e9
+        print(" [az=%5.1f° el=%5.1f° dist=%5dkm/%4.1fms ds~%6dHz]"%(az.degrees,el.degrees,dist.km,dist.m/ 299792.458,shift),end="")
+        if args.verbose:
+            print(" [ar=%+5.2f°/s er=%+5.2f°/s dr=%+5.2fkm/s]"%(az_r.degrees.per_second,el_r.degrees.per_second,dist_r.km_per_s),end="")
         print("")
         for t,o in out:
             print(t.utc_datetime().astimezone().strftime("%H:%M:%S"),"-",deltatstr(t.utc_datetime()), o)
 
-satev=generate_events(satlist)
-nice=format_events(satev, nitems=items)
 
-try:
-    fd=None
-    old_term=None
-    if repeat:
+if __name__ == "__main__":
+    parse_args()
+    observer=read_observer(args.loc)
+    recpos=Topos(latitude_degrees=observer['lat'], longitude_degrees=observer['lon'], elevation_m=observer['alt'])
+    satlist = loadTLE(args.tlefile)
+    now=datetime.now(timezone.utc)
+    tnow = mytsutc(now)
+
+    sat0=satlist[0]
+    days = tnow - sat0.epoch
+    if args.verbose:
+        print('TLE file is %.2f days old'%days)
+
+    if abs(days)>5:
+        print('SAT EPOCH too far away. (%.2f days)'%days, file=sys.stderr)
+        sys.exit(-1)
+
+    satev=generate_events(satlist)
+    nice=format_events(satev, nitems=args.items)
+
+    if not args.repeat:
+        print_events(nice)
+        exit(0)
+
+    try:
+        fd=None
+        old_term=None
         import curses
         curses.initscr()
         curses.endwin()
@@ -199,17 +262,16 @@ def print_events(nice):
         new_term[3] = (new_term[3] & ~termios.ICANON & ~termios.ECHO)
         termios.tcsetattr(fd, termios.TCSAFLUSH, new_term)
 
-    while repeat:
-        now=datetime.now(timezone.utc)
-        if repeat:
+        while True:
+            items=args.items
+            now=datetime.now(timezone.utc)
             sys.stdout.buffer.write(curses.tigetstr("home"))
-            print(os.path.basename(__file__)," - ",now.astimezone().strftime("%H:%M:%S"))
-        print_events(nice)
-        if nice[-1][1][-1][0].utc_datetime()<now: # Exit when last sat in list has set.
-            print()
-            print("done.")
-            break
-        if repeat:
+            print(os.path.basename(__file__)," - ",now.astimezone().strftime("%H:%M:%S")," - ",args.loc)
+            print_events(nice)
+            if nice[-1][1][-1][0].utc_datetime()<now: # Exit when last sat in list has set.
+                print()
+                print("done.")
+                break
 #            sys.stdout.buffer.write(curses.tigetstr("el"))
             sys.stdout.flush()
             # sleep a bit if no input
@@ -231,10 +293,10 @@ def print_events(nice):
                     print("Unknown key pressed")
                     select.select([sys.stdin], [], [], 1)
 
-except KeyboardInterrupt:
-    pass
+    except KeyboardInterrupt:
+        print("^C")
+        pass
 
-finally:
-    if repeat:
+    finally:
         sys.stdout.buffer.write(curses.tigetstr("cnorm"))
         termios.tcsetattr(fd, termios.TCSAFLUSH, old_term)