From 550399a3605a1fea45978194513af5bd51a4a2bc Mon Sep 17 00:00:00 2001 From: Pierros Papadeas Date: Sun, 12 May 2019 17:11:15 +0300 Subject: [PATCH] Add calibrate.py for astrometric solving --- acquire.py | 3 +- calibrate.py | 51 +++++++++ configuration.ini-dist | 3 + imgstat.py | 51 +++++---- keogram.py | 31 ++++-- process.py | 54 +++++----- setup.cfg | 4 + stvid/astrometry.py | 1 + stvid/extract.py | 233 +++++++++++++++++++---------------------- stvid/satellite.py | 35 ++++--- stvid/stio.py | 117 ++++++++++----------- stvid/utils.py | 92 ++++++++-------- update_tle.py | 1 - 13 files changed, 367 insertions(+), 309 deletions(-) create mode 100755 calibrate.py create mode 100644 setup.cfg diff --git a/acquire.py b/acquire.py index 773deac..8fdc73e 100755 --- a/acquire.py +++ b/acquire.py @@ -205,7 +205,6 @@ if __name__ == '__main__': logging.basicConfig(filename=os.path.join(path, "acquire.log"), level=logging.DEBUG) - # Set location loc = EarthLocation(lat=cfg.getfloat('Common', 'observer_lat')*u.deg, lon=cfg.getfloat('Common', 'observer_lon')*u.deg, @@ -218,7 +217,7 @@ if __name__ == '__main__': # Get sunrise and sunset times state, tset, trise = get_sunset_and_sunrise(tnow, loc, refalt_set, refalt_rise) - + # Start/end logic if state == "sun never rises": logging.info("The sun never rises. Exiting program.") diff --git a/calibrate.py b/calibrate.py new file mode 100755 index 0000000..6b09245 --- /dev/null +++ b/calibrate.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python +from __future__ import print_function +import os +import configparser +import argparse +import subprocess + +if __name__ == '__main__': + # Read commandline options + conf_parser = argparse.ArgumentParser(description='Plate solve FITS file' + + ' and add WCS on header') + conf_parser.add_argument("-c", "--conf_file", + help="Specify configuration file. If no file" + + " is specified 'configuration.ini' is used.", + metavar="FILE") + conf_parser.add_argument("-d", "--directory", + help="Specify directory of observations. If no" + + " directory is specified parent will be used.", + metavar='DIR', dest='file_dir', default=".") + + args = conf_parser.parse_args() + + # Process commandline options and parse configuration + cfg = configparser.ConfigParser(inline_comment_prefixes=('#', ';')) + if args.conf_file: + cfg.read([args.conf_file]) + else: + cfg.read('configuration.ini') + + path = args.file_dir + extension = 'fits' + files = sorted((f for f in os.listdir(path) if f.endswith('.' + extension)), + key=lambda x: os.path.getctime(os.path.join(path, x))) + + if files: + file_for_astrometry = os.path.join(path, files[0]) + print("Found " + file_for_astrometry + " for astrometric solving.") + + sex_config = cfg.get('Astrometry', 'sex_config') + + # Format command + command = "solve-field -O -y -u app -L 37 -H 40 --use-sextractor " + \ + "--sextractor-config %s --downsample 2 --x-column X_IMAGE " % sex_config + \ + "--y-column Y_IMAGE --sort-column MAG_AUTO --sort-ascending " + \ + "%s" % file_for_astrometry + + # Run sextractor + subprocess.run(command, shell=True, stderr=subprocess.STDOUT) + + else: + print("No fits file found for astrometric solving.") diff --git a/configuration.ini-dist b/configuration.ini-dist index 6b67de5..67f8a95 100644 --- a/configuration.ini-dist +++ b/configuration.ini-dist @@ -24,3 +24,6 @@ device_id = 0 camera_x = 720 # Camera horizontal pixel count camera_y = 576 # Camera vertical pixel count camera_frames = 250 # Camera frames for each image (25fps*10sec=250) + +[Astrometry] +sex_config = /path/to/solve.sex diff --git a/imgstat.py b/imgstat.py index 95170f8..81ca2e2 100755 --- a/imgstat.py +++ b/imgstat.py @@ -1,6 +1,5 @@ #!/usr/bin/env python3 from __future__ import print_function -import numpy as np from astropy.io import ascii import matplotlib.pyplot as plt import matplotlib.dates as mdates @@ -14,34 +13,43 @@ import os if __name__ == "__main__": # Read commandline options conf_parser = argparse.ArgumentParser(description='Plot image statistics') - conf_parser.add_argument("-c", "--conf_file", + conf_parser.add_argument("-c", + "--conf_file", help="Specify configuration file. If no file" + " is specified 'configuration.ini' is used.", metavar="FILE") - conf_parser.add_argument("-i", "--input", + conf_parser.add_argument("-i", + "--input", help="Specify file to be processed. If no file" + " is specified ./imgstat.csv will be used.", - metavar='FILE', default="./imgstat.csv") - conf_parser.add_argument("-d", "--directory", + metavar='FILE', + default="./imgstat.csv") + conf_parser.add_argument("-d", + "--directory", help="Specify directory of observations. If no" + " directory is specified parent will be used.", - metavar='DIR', dest='file_dir', default=".") - conf_parser.add_argument("-o", "--output", - help="Specify output file. Default is 'imgstat.png'", - metavar='FILE', default="./imgstat.png") - + metavar='DIR', + dest='file_dir', + default=".") + conf_parser.add_argument( + "-o", + "--output", + help="Specify output file. Default is 'imgstat.png'", + metavar='FILE', + default="./imgstat.png") + args = conf_parser.parse_args() - + # Process commandline options and parse configuration cfg = configparser.ConfigParser(inline_comment_prefixes=('#', ';')) if args.conf_file: cfg.read([args.conf_file]) else: cfg.read('configuration.ini') - + # Move to processing directory os.chdir(args.file_dir) - + table = ascii.read(args.input, format="csv") t = Time(table['mjd'], format="mjd", scale="utc") @@ -49,13 +57,16 @@ if __name__ == "__main__": pos = SkyCoord(ra=table['ra'], dec=table['de'], frame="icrs", unit="deg") # Set location - loc = EarthLocation(lat=cfg.getfloat('Common', 'observer_lat')*u.deg, - lon=cfg.getfloat('Common', 'observer_lon')*u.deg, - height=cfg.getfloat('Common', 'observer_el')*u.m) - + loc = EarthLocation(lat=cfg.getfloat('Common', 'observer_lat') * u.deg, + lon=cfg.getfloat('Common', 'observer_lon') * u.deg, + height=cfg.getfloat('Common', 'observer_el') * u.m) + pa = pos.transform_to(AltAz(obstime=t, location=loc)) - fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(20, 10), sharex=True) + fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, + 1, + figsize=(20, 10), + sharex=True) date_format = mdates.DateFormatter("%F\n%H:%M:%S") fig.autofmt_xdate(rotation=0, ha="center") @@ -85,8 +96,6 @@ if __name__ == "__main__": ax4.set_xlabel("Time (UTC)") ax4.set_ylabel("Residuals (arcseconds)") ax4.legend() - + plt.tight_layout() plt.savefig(args.output) - - diff --git a/keogram.py b/keogram.py index 10ff9b0..02013c5 100755 --- a/keogram.py +++ b/keogram.py @@ -10,6 +10,7 @@ import matplotlib.pyplot as plt import matplotlib.dates as mdates from astropy.time import Time + def generate_keogram(path): # Get files fnames = sorted(glob.glob(os.path.join(path, "processed/2*.fits"))) @@ -17,13 +18,13 @@ def generate_keogram(path): # Allocate arrays nx = len(fnames) ny = cfg.getint('Camera', 'camera_y') - ixmid = cfg.getint('Camera', 'camera_x')//2 - keogram = np.zeros(nx*ny).reshape(ny, nx) + ixmid = cfg.getint('Camera', 'camera_x') // 2 + keogram = np.zeros(nx * ny).reshape(ny, nx) mjds = np.zeros(nx) - + # Get data for i, fname in enumerate(fnames): - if i%10==0: + if i % 10 == 0: print(i, fname) # Read file @@ -41,15 +42,19 @@ def generate_keogram(path): if __name__ == "__main__": # Read commandline options conf_parser = argparse.ArgumentParser(description='Process captured' + - ' video frames.') - conf_parser.add_argument("-c", "--conf_file", + ' video frames.') + conf_parser.add_argument("-c", + "--conf_file", help="Specify configuration file. If no file" + " is specified 'configuration.ini' is used.", metavar="FILE") - conf_parser.add_argument("-d", "--directory", + conf_parser.add_argument("-d", + "--directory", help="Specify directory of observations. If no" + " directory is specified parent will be used.", - metavar='DIR', dest='file_dir', default=".") + metavar='DIR', + dest='file_dir', + default=".") args = conf_parser.parse_args() @@ -71,16 +76,20 @@ if __name__ == "__main__": # Time limits tmin = mdates.date2num(Time(np.min(mjds), format="mjd").datetime) tmax = mdates.date2num(Time(np.max(mjds), format="mjd").datetime) - + # Plot keogram fig, ax = plt.subplots(figsize=(15, 5)) - ax.imshow(np.log10(keogram), aspect="auto", origin="lower", cmap="magma", extent=[tmin, tmax, 0, 1]) + ax.imshow(np.log10(keogram), + aspect="auto", + origin="lower", + cmap="magma", + extent=[tmin, tmax, 0, 1]) ax.axes.get_yaxis().set_visible(False) ax.xaxis_date() date_format = mdates.DateFormatter("%F\n%H:%M:%S") ax.xaxis.set_major_formatter(date_format) fig.autofmt_xdate(rotation=0, ha="center") - + plt.tight_layout() plt.savefig(os.path.join(args.file_dir, "keogram.png")) diff --git a/process.py b/process.py index 6df489a..c18164e 100755 --- a/process.py +++ b/process.py @@ -23,15 +23,19 @@ import shutil if __name__ == "__main__": # Read commandline options conf_parser = argparse.ArgumentParser(description='Process captured' + - ' video frames.') - conf_parser.add_argument("-c", "--conf_file", + ' video frames.') + conf_parser.add_argument("-c", + "--conf_file", help="Specify configuration file. If no file" + " is specified 'configuration.ini' is used.", metavar="FILE") - conf_parser.add_argument("-d", "--directory", + conf_parser.add_argument("-d", + "--directory", help="Specify directory of observations. If no" + " directory is specified parent will be used.", - metavar='DIR', dest='file_dir', default=".") + metavar='DIR', + dest='file_dir', + default=".") args = conf_parser.parse_args() @@ -47,9 +51,9 @@ if __name__ == "__main__": warnings.simplefilter("ignore", AstropyWarning) # Set location - loc = EarthLocation(lat=cfg.getfloat('Common', 'observer_lat')*u.deg, - lon=cfg.getfloat('Common', 'observer_lon')*u.deg, - height=cfg.getfloat('Common', 'observer_el')*u.m) + loc = EarthLocation(lat=cfg.getfloat('Common', 'observer_lat') * u.deg, + lon=cfg.getfloat('Common', 'observer_lon') * u.deg, + height=cfg.getfloat('Common', 'observer_el') * u.m) # Extract settings # Minimum predicted velocity (pixels/s) @@ -70,7 +74,7 @@ if __name__ == "__main__": # Statistics file fstat = open("imgstat.csv", "w") fstat.write("fname,mjd,ra,de,rmsx,rmsy,mean,std,nstars,nused\n") - + # Create output dirs path = args.file_dir results_path = os.path.join(cfg.get('Common', 'results_path'), @@ -88,53 +92,55 @@ if __name__ == "__main__": while True: # Get files fnames = sorted(glob.glob("2*.fits")) - + # Loop over files for fname in fnames: # Generate star catalog pix_catalog = generate_star_catalog(fname) # Calibrate astrometry - calibrate_from_reference(fname, "test.fits", - pix_catalog) + calibrate_from_reference(fname, "test.fits", pix_catalog) # Generate satellite predictions generate_satellite_predictions(fname) - + # Detect lines with 3D Hough transform ids = find_hough3d_lines(fname) # Extract tracks - extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, results_path) - + extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, + results_path) + # Stars available and used nused = np.sum(pix_catalog.flag == 1) nstars = pix_catalog.nstars - + # Get properties ff = fourframe(fname) # Write output - output = "%s %10.6f %10.6f %4d/%4d %5.1f %5.1f %6.2f +- %6.2f"%(ff.fname, ff.crval[0], ff.crval[1], nused, nstars, 3600.0*ff.crres[0], 3600.0*ff.crres[1], np.mean(ff.zavg), np.std(ff.zavg)) + output = "%s %10.6f %10.6f %4d/%4d %5.1f %5.1f %6.2f +- %6.2f" % ( + ff.fname, ff.crval[0], ff.crval[1], nused, nstars, + 3600.0 * ff.crres[0], 3600.0 * ff.crres[1], np.mean( + ff.zavg), np.std(ff.zavg)) if is_calibrated(ff): print(colored(output, "green")) else: print(colored(output, "red")) - fstat.write(("%s,%.8lf,%.6f,%.6f,%.3f,%.3f,%.3f," + - "%.3f,%d,%d\n") % (ff.fname, ff.mjd, ff.crval[0], - ff.crval[1], 3600*ff.crres[0], - 3600*ff.crres[1], np.mean(ff.zavg), - np.std(ff.zavg), nstars, nused)) + fstat.write( + ("%s,%.8lf,%.6f,%.6f,%.3f,%.3f,%.3f," + "%.3f,%d,%d\n") % + (ff.fname, ff.mjd, ff.crval[0], ff.crval[1], + 3600 * ff.crres[0], 3600 * ff.crres[1], np.mean( + ff.zavg), np.std(ff.zavg), nstars, nused)) # Move processed files shutil.move(fname, "processed") shutil.move(fname + ".png", "processed") shutil.move(fname + ".id", "processed") shutil.move(fname + ".cat", "processed") - + # Sleep time.sleep(10) - + # Close files fstat.close() - diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..0b0122b --- /dev/null +++ b/setup.cfg @@ -0,0 +1,4 @@ +[flake8] +max-complexity = 25 +max-line-length = 99 +ignore = F403, W504, E226, W503 diff --git a/stvid/astrometry.py b/stvid/astrometry.py index 28e6324..03875f2 100644 --- a/stvid/astrometry.py +++ b/stvid/astrometry.py @@ -207,6 +207,7 @@ def calibrate_from_reference(fname, ref, pix_catalog): return w, rmsx, rmsy + def is_calibrated(ff): if (3600.0*ff.crres[0] < 1e-3) | \ (3600.0*ff.crres[1] < 1e-3) | \ diff --git a/stvid/extract.py b/stvid/extract.py index ada15a5..c16c3da 100644 --- a/stvid/extract.py +++ b/stvid/extract.py @@ -1,7 +1,6 @@ #!/usr/bin/env python from __future__ import print_function import os -import glob import shutil from stvid.stio import fourframe, satid, observation from stvid.astrometry import is_calibrated @@ -11,33 +10,36 @@ from scipy import optimize, ndimage from termcolor import colored import datetime + # Gaussian model def model(a, nx, ny): x, y = np.meshgrid(np.arange(nx), np.arange(ny)) - dx, dy = (x-a[0])/a[2], (y-a[1])/a[2] - arg = -0.5*(dx**2+dy**2) - return a[3]*np.exp(arg)+a[4] + dx, dy = (x - a[0]) / a[2], (y - a[1]) / a[2] + arg = -0.5 * (dx**2 + dy**2) + return a[3] * np.exp(arg) + a[4] + # Residual function def residual(a, img): ny, nx = img.shape mod = model(a, nx, ny) - return (img-mod).ravel() + return (img - mod).ravel() + # Find peak def peakfind(img, w=1.0): # Find approximate location ny, nx = img.shape i = np.argmax(img) - y0 = int(i/nx) - x0 = i-y0*nx + y0 = int(i / nx) + x0 = i - y0 * nx # Image properties imgavg = np.mean(img) imgstd = np.std(img) # Estimate - a = np.array([x0, y0, w, img[y0, x0]-imgavg, imgavg]) + a = np.array([x0, y0, w, img[y0, x0] - imgavg, imgavg]) q, cov_q, infodict, mesg, ier = optimize.leastsq(residual, a, args=(img), @@ -47,82 +49,87 @@ def peakfind(img, w=1.0): xc, yc, w = q[0], q[1], q[2] # Significance - sigma = (a[3]-imgavg)/(imgstd+1e-9) + sigma = (a[3] - imgavg) / (imgstd + 1e-9) return xc, yc, w, sigma + # Plot selection def plot_selection(id, x0, y0, dt=2.0, w=10.0): - dx, dy = id.x1-id.x0, id.y1-id.y0 + dx, dy = id.x1 - id.x0, id.y1 - id.y0 ang = np.arctan2(dy, dx) - r = np.sqrt(dx**2+dy**2) - drdt = r/(id.t1-id.t0) + r = np.sqrt(dx**2 + dy**2) + drdt = r / (id.t1 - id.t0) sa, ca = np.sin(ang), np.cos(ang) - dx = np.array([-dt, -dt, dt, dt, -dt])*drdt + dx = np.array([-dt, -dt, dt, dt, -dt]) * drdt dy = np.array([w, -w, -w, w, w]) - x = ca*dx-sa*dy+x0 - y = sa*dx+ca*dy+y0 + x = ca * dx - sa * dy + x0 + y = sa * dx + ca * dy + y0 ppg.pgsci(7) ppg.pgline(x, y) return + # Check if point is inside selection def inside_selection(ident, tmid, x0, y0, dt=2.0, w=10.0): - dx, dy = ident.x1-ident.x0, ident.y1-ident.y0 + dx, dy = ident.x1 - ident.x0, ident.y1 - ident.y0 ang = -np.arctan2(dy, dx) - r = np.sqrt(dx**2+dy**2) - drdt = r/(ident.t1-ident.t0) + r = np.sqrt(dx**2 + dy**2) + drdt = r / (ident.t1 - ident.t0) sa, ca = np.sin(ang), np.cos(ang) - xmid = ident.x0+ident.dxdt*tmid - ymid = ident.y0+ident.dydt*tmid - - dx, dy = x0-xmid, y0-ymid - rm = ca*dx-sa*dy - wm = sa*dx+ca*dy - dtm = rm/drdt + xmid = ident.x0 + ident.dxdt * tmid + ymid = ident.y0 + ident.dydt * tmid + + dx, dy = x0 - xmid, y0 - ymid + rm = ca * dx - sa * dy + wm = sa * dx + ca * dy + dtm = rm / drdt if (abs(wm) < w) & (abs(dtm) < dt): return True else: return False + # Get COSPAR ID def get_cospar(norad, nfd): f = open(os.path.join(os.getenv("ST_DATADIR"), "data/desig.txt")) lines = f.readlines() f.close() - + try: cospar = ([line for line in lines if str(norad) in line])[0].split()[1] except IndexError: t = datetime.datetime.strptime(nfd[:-4], "%Y-%m-%dT%H:%M:%S") - doy = int(t.strftime("%y%j"))+500 - cospar = "%sA"%doy + doy = int(t.strftime("%y%j")) + 500 + cospar = "%sA" % doy return "%2s %s" % (cospar[0:2], cospar[2:]) + # IOD position format 2: RA/DEC = HHMMmmm+DDMMmm MX (MX in minutes of arc) def format_position(ra, de): - ram = 60.0*ra/15.0 - rah = int(np.floor(ram/60.0)) - ram -= 60.0*rah + ram = 60.0 * ra / 15.0 + rah = int(np.floor(ram / 60.0)) + ram -= 60.0 * rah des = np.sign(de) - dem = 60.0*np.abs(de) - ded = int(np.floor(dem/60.0)) - dem -= 60.0*ded + dem = 60.0 * np.abs(de) + ded = int(np.floor(dem / 60.0)) + dem -= 60.0 * ded if des == -1: sign = "-" else: sign = "+" - return ("%02d%06.3f%c%02d%05.2f" - % (rah, ram, sign, ded, dem)).replace(".", "") + return ("%02d%06.3f%c%02d%05.2f" % (rah, ram, sign, ded, dem)).replace( + ".", "") + # Format IOD line def format_iod_line(norad, cospar, site_id, t, ra, de): @@ -132,12 +139,10 @@ def format_iod_line(norad, cospar, site_id, t, ra, de): .replace(":", "") \ .replace(".", "") - return "%05d %-9s %04d G %s 17 25 %s 37 S" % (norad, - cospar, - site_id, - tstr, + return "%05d %-9s %04d G %s 17 25 %s 37 S" % (norad, cospar, site_id, tstr, pstr) + def store_results(ident, fname, path, iod_line): # Find destination if ident.catalog.find("classfd.tle") > 0: @@ -159,31 +164,32 @@ def store_results(ident, fname, path, iod_line): # Print iod line print(colored(iod_line, color)) - + # Copy files shutil.copy2(fname, dest) shutil.copy2(fname + ".cat", dest) shutil.copy2(fname + ".id", dest) shutil.copy2(fname + ".png", dest) try: - shutil.move(fname.replace(".fits", "_%05d.png"%ident.norad), dest) - except: + shutil.move(fname.replace(".fits", "_%05d.png" % ident.norad), dest) + except Exception: pass - + # Write iodline fp = open(outfname, "a") fp.write("%s\n" % iod_line) fp.close() - + return + def plot_header(fname, ff, iod_line): # ppgplot arrays heat_l = np.array([0.0, 0.2, 0.4, 0.6, 1.0]) heat_r = np.array([0.0, 0.5, 1.0, 1.0, 1.0]) heat_g = np.array([0.0, 0.0, 0.5, 1.0, 1.0]) heat_b = np.array([0.0, 0.0, 0.0, 0.3, 1.0]) - + # Plot ppg.pgopen(fname) ppg.pgpap(0.0, 1.0) @@ -191,29 +197,21 @@ def plot_header(fname, ff, iod_line): ppg.pgsch(0.8) ppg.pgmtxt("T", 6.0, 0.0, 0.0, - "UT Date: %.23s COSPAR ID: %04d" - % (ff.nfd, ff.site_id)) + "UT Date: %.23s COSPAR ID: %04d" % (ff.nfd, ff.site_id)) if is_calibrated(ff): ppg.pgsci(1) else: ppg.pgsci(2) - ppg.pgmtxt("T", 4.8, 0.0, 0.0, - "R.A.: %10.5f (%4.1f'') Decl.: %10.5f (%4.1f'')" - % (ff.crval[0], - 3600.0*ff.crres[0], - ff.crval[1], - 3600.0*ff.crres[1])) + ppg.pgmtxt( + "T", 4.8, 0.0, 0.0, "R.A.: %10.5f (%4.1f'') Decl.: %10.5f (%4.1f'')" % + (ff.crval[0], 3600.0 * ff.crres[0], ff.crval[1], 3600.0 * ff.crres[1])) ppg.pgsci(1) - ppg.pgmtxt("T", 3.6, 0.0, 0.0, - ("FoV: %.2f\\(2218)x%.2f\\(2218) " - "Scale: %.2f''x%.2f'' pix\\u-1\\d") - % (ff.wx, ff.wy, 3600.0*ff.sx, 3600.0*ff.sy)) - ppg.pgmtxt("T", 2.4, 0.0, 0.0, - "Stat: %5.1f+-%.1f (%.1f-%.1f)" - % (np.mean(ff.zmax), - np.std(ff.zmax), - ff.zmaxmin, - ff.zmaxmax)) + ppg.pgmtxt("T", 3.6, 0.0, 0.0, ("FoV: %.2f\\(2218)x%.2f\\(2218) " + "Scale: %.2f''x%.2f'' pix\\u-1\\d") % + (ff.wx, ff.wy, 3600.0 * ff.sx, 3600.0 * ff.sy)) + ppg.pgmtxt( + "T", 2.4, 0.0, 0.0, "Stat: %5.1f+-%.1f (%.1f-%.1f)" % + (np.mean(ff.zmax), np.std(ff.zmax), ff.zmaxmin, ff.zmaxmax)) ppg.pgmtxt("T", 0.3, 0.0, 0.0, iod_line) ppg.pgsch(1.0) @@ -228,14 +226,14 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): ff = fourframe(fname) # Skip saturated frames - if np.sum(ff.zavg > 240.0)/float(ff.nx*ff.ny) > 0.95: + if np.sum(ff.zavg > 240.0) / float(ff.nx * ff.ny) > 0.95: return # Read satelite IDs try: - f = open(fname+".id") + f = open(fname + ".id") except OSError: - print("Cannot open", fname+".id") + print("Cannot open", fname + ".id") else: lines = f.readlines() f.close() @@ -244,7 +242,7 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): # Parse identifications idents = [satid(line) for line in lines] - + # Identify unknowns for ident0 in idents: if ident0.catalog == "unidentified": @@ -264,67 +262,60 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): ident0.state = ident1.state ident1.state = "remove" break - + # Loop over identifications for ident in idents: # Skip superseded unknowns if ident.state == "remove": continue - + # Skip slow moving objects - drdt = np.sqrt(ident.dxdt**2+ident.dydt**2) + drdt = np.sqrt(ident.dxdt**2 + ident.dydt**2) if drdt < drdtmin: continue # Extract significant pixels along a track - x, y, t, sig = ff.significant_pixels_along_track(trksig, - ident.x0, - ident.y0, - ident.dxdt, - ident.dydt, - trkrmin) + x, y, t, sig = ff.significant_pixels_along_track( + trksig, ident.x0, ident.y0, ident.dxdt, ident.dydt, trkrmin) # Fit tracks if len(t) > ntrkmin: # Get times tmin = np.min(t) tmax = np.max(t) - tmid = 0.5*(tmax+tmin) - mjd = ff.mjd+tmid/86400.0 + tmid = 0.5 * (tmax + tmin) + mjd = ff.mjd + tmid / 86400.0 # Skip if no variance in time - if np.std(t-tmid) == 0.0: + if np.std(t - tmid) == 0.0: continue # Very simple polynomial fit; no weighting, no cleaning - px = np.polyfit(t-tmid, x, 1) - py = np.polyfit(t-tmid, y, 1) + px = np.polyfit(t - tmid, x, 1) + py = np.polyfit(t - tmid, y, 1) # Extract results x0, y0 = px[1], py[1] dxdt, dydt = px[0], py[0] - xmin = x0+dxdt*(tmin-tmid) - ymin = y0+dydt*(tmin-tmid) - xmax = x0+dxdt*(tmax-tmid) - ymax = y0+dydt*(tmax-tmid) + xmin = x0 + dxdt * (tmin - tmid) + ymin = y0 + dydt * (tmin - tmid) + xmax = x0 + dxdt * (tmax - tmid) + ymax = y0 + dydt * (tmax - tmid) cospar = get_cospar(ident.norad, ff.nfd) obs = observation(ff, mjd, x0, y0) - iod_line = "%s" % format_iod_line(ident.norad, - cospar, - ff.site_id, - obs.nfd, - obs.ra, - obs.de) + iod_line = "%s" % format_iod_line(ident.norad, cospar, ff.site_id, + obs.nfd, obs.ra, obs.de) # Create diagnostic plot - plot_header(fname.replace(".fits", "_%05d.png/png" % ident.norad), ff, iod_line) + plot_header(fname.replace(".fits", "_%05d.png/png" % ident.norad), + ff, iod_line) - ppg.pgimag(ff.zmax, ff.nx, ff.ny, 0, ff.nx-1, - 0, ff.ny-1, ff.zmaxmax, ff.zmaxmin, tr) + ppg.pgimag(ff.zmax, ff.nx, ff.ny, 0, ff.nx - 1, 0, ff.ny - 1, + ff.zmaxmax, ff.zmaxmin, tr) ppg.pgbox("BCTSNI", 0., 0, "BCTSNI", 0., 0) ppg.pgstbg(1) - + ppg.pgsci(0) if ident.catalog.find("classfd.tle") > 0: ppg.pgsci(4) @@ -338,7 +329,7 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): ppg.pgtext(np.array([x0]), np.array([y0]), " %05d" % ident.norad) ppg.pgsch(1.0) ppg.pgsci(1) - + ppg.pgend() # Store results @@ -347,7 +338,7 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): elif ident.catalog.find("classfd.tle") > 0: # Track and stack t = np.linspace(0.0, ff.texp) - x, y = ident.x0+ident.dxdt*t, ident.y0+ident.dydt*t + x, y = ident.x0 + ident.dxdt * t, ident.y0 + ident.dydt * t c = (x > 0) & (x < ff.nx) & (y > 0) & (y < ff.ny) # Skip if no points selected @@ -356,27 +347,27 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): # Compute track tmid = np.mean(t[c]) - mjd = ff.mjd+tmid/86400.0 - xmid = ident.x0+ident.dxdt*tmid - ymid = ident.y0+ident.dydt*tmid - ztrk = ndimage.gaussian_filter(ff.track(ident.dxdt, ident.dydt, tmid), - 1.0) - vmin = np.mean(ztrk)-2.0*np.std(ztrk) - vmax = np.mean(ztrk)+6.0*np.std(ztrk) + mjd = ff.mjd + tmid / 86400.0 + xmid = ident.x0 + ident.dxdt * tmid + ymid = ident.y0 + ident.dydt * tmid + ztrk = ndimage.gaussian_filter( + ff.track(ident.dxdt, ident.dydt, tmid), 1.0) + vmin = np.mean(ztrk) - 2.0 * np.std(ztrk) + vmax = np.mean(ztrk) + 6.0 * np.std(ztrk) # Select region - xmin = int(xmid-100) - xmax = int(xmid+100) - ymin = int(ymid-100) - ymax = int(ymid+100) + xmin = int(xmid - 100) + xmax = int(xmid + 100) + ymin = int(ymid - 100) + ymax = int(ymid + 100) if xmin < 0: xmin = 0 if ymin < 0: ymin = 0 if xmax > ff.nx: - xmax = ff.nx-1 + xmax = ff.nx - 1 if ymax > ff.ny: - ymax = ff.ny-1 + ymax = ff.ny - 1 # Find peak x0, y0, w, sigma = peakfind(ztrk[ymin:ymax, xmin:xmax]) @@ -394,20 +385,15 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): # Format IOD line cospar = get_cospar(ident.norad, ff.nfd) obs = observation(ff, mjd, x0, y0) - iod_line = "%s" % format_iod_line(ident.norad, - cospar, - ff.site_id, - obs.nfd, - obs.ra, - obs.de) - + iod_line = "%s" % format_iod_line(ident.norad, cospar, ff.site_id, + obs.nfd, obs.ra, obs.de) # Create diagnostic plot - plot_header(fname.replace(".fits", "_%05d.png/png" % ident.norad), ff, iod_line) + plot_header(fname.replace(".fits", "_%05d.png/png" % ident.norad), + ff, iod_line) - - ppg.pgimag(ztrk, ff.nx, ff.ny, 0, ff.nx-1, - 0, ff.ny-1, vmax, vmin, tr) + ppg.pgimag(ztrk, ff.nx, ff.ny, 0, ff.nx - 1, 0, ff.ny - 1, vmax, + vmin, tr) ppg.pgbox("BCTSNI", 0., 0, "BCTSNI", 0., 0) ppg.pgstbg(1) @@ -423,8 +409,7 @@ def extract_tracks(fname, trkrmin, drdtmin, trksig, ntrkmin, path): ppg.pgdraw(ident.x1, ident.y1) ppg.pgpt(np.array([x0]), np.array([y0]), 4) ppg.pgsch(0.65) - ppg.pgtext(np.array([ident.x0]), - np.array([ident.y0]), + ppg.pgtext(np.array([ident.x0]), np.array([ident.y0]), " %05d" % ident.norad) ppg.pgsch(1.0) ppg.pgsci(1) diff --git a/stvid/satellite.py b/stvid/satellite.py index 740775b..bcc1c22 100644 --- a/stvid/satellite.py +++ b/stvid/satellite.py @@ -1,17 +1,17 @@ #!/usr/bin/env python from __future__ import print_function -import os import subprocess from stvid.stio import fourframe from stvid.stio import satid -from stvid.stio import observation + def generate_satellite_predictions(fname): # Format command command = "satid %s %s.png/png" % (fname, fname) # Run command - output = subprocess.check_output(command, shell=True, + output = subprocess.check_output(command, + shell=True, stderr=subprocess.STDOUT) return output @@ -28,9 +28,9 @@ def find_hough3d_lines(fname, ntrkmin=20, dr=8): x, y, z, t, sig = ff.selection_mask(5.0, 40.0) # Skip if not enough points - if len(t)<2: + if len(t) < 2: return [] - + # Save points to temporary file with open("/tmp/hough.dat", "w") as f: for i in range(len(t)): @@ -38,12 +38,13 @@ def find_hough3d_lines(fname, ntrkmin=20, dr=8): f.close() # Run 3D Hough line-finding algorithm - command = "hough3dlines -dx %d -minvotes %d %s" % (dr, ntrkmin, "/tmp/hough.dat") + command = "hough3dlines -dx %d -minvotes %d %s" % (dr, ntrkmin, + "/tmp/hough.dat") try: output = subprocess.check_output(command, shell=True, stderr=subprocess.STDOUT) - except: + except Exception: return [] # Clean output (a bit cluncky) @@ -56,16 +57,18 @@ def find_hough3d_lines(fname, ntrkmin=20, dr=8): # Generate identifications lines = [] s = cleaned_output.split() - for i in range(len(s)//7): + for i in range(len(s) // 7): # Extract points - x0, y0, z0 = float(s[1+7*i]), float(s[2+7*i]), float(s[3+7*i]) - dx, dy, dz = float(s[4+7*i]), float(s[5+7*i]), float(s[6+7*i]) - + x0, y0, z0 = float(s[1 + 7 * i]), float(s[2 + 7 * i]), float(s[3 + + 7 * i]) + dx, dy, dz = float(s[4 + 7 * i]), float(s[5 + 7 * i]), float(s[6 + + 7 * i]) + # Reconstruct start and end points - xmin = x0-z0*dx/(dz+1e-9) - xmax = x0+(ff.nz-z0)*dx/(dz+1e-9) - ymin = y0-z0*dy/(dz+1e-9) - ymax = y0+(ff.nz-z0)*dy/(dz+1e-9) + xmin = x0 - z0 * dx / (dz + 1e-9) + xmax = x0 + (ff.nz - z0) * dx / (dz + 1e-9) + ymin = y0 - z0 * dy / (dz + 1e-9) + ymax = y0 + (ff.nz - z0) * dy / (dz + 1e-9) # Output line line = "%.23s %8.3f %8.3f %8.3f %8.3f %8.5f %s unidentified sunlit\n" %\ @@ -77,5 +80,5 @@ def find_hough3d_lines(fname, ntrkmin=20, dr=8): for line in lines: fp.write(line) fp.close() - + return [satid(line) for line in lines] diff --git a/stvid/stio.py b/stvid/stio.py index f419422..8cfb014 100644 --- a/stvid/stio.py +++ b/stvid/stio.py @@ -22,7 +22,9 @@ class observation: self.nfd = Time(self.mjd, format='mjd', scale='utc').isot # Correct for rotation - tobs = Time(ff.mjd+0.5*ff.texp/86400.0, format='mjd', scale='utc') + tobs = Time(ff.mjd + 0.5 * ff.texp / 86400.0, + format='mjd', + scale='utc') tobs.delta_ut1_utc = 0 hobs = tobs.sidereal_time("mean", longitude=0.0).degree tmid = Time(self.mjd, format='mjd', scale='utc') @@ -31,7 +33,7 @@ class observation: # Compute ra/dec world = ff.w.wcs_pix2world(np.array([[self.x0, self.y0]]), 1) - self.ra = world[0, 0]+hobs-hmid + self.ra = world[0, 0] + hobs - hmid self.de = world[0, 1] @@ -50,15 +52,13 @@ class satid: self.norad = int(s[6]) self.catalog = s[7] self.state = s[8] - self.dxdt = (self.x1-self.x0)/(self.t1-self.t0) - self.dydt = (self.y1-self.y0)/(self.t1-self.t0) + self.dxdt = (self.x1 - self.x0) / (self.t1 - self.t0) + self.dydt = (self.y1 - self.y0) / (self.t1 - self.t0) def __repr__(self): - return "%s %f %f %f -> %f %f %f %d %s %s" % (self.nfd, self.x0, - self.y0, self.t0, - self.x1, self.y1, - self.t1, self.norad, - self.catalog, self.state) + return "%s %f %f %f -> %f %f %f %d %s %s" % ( + self.nfd, self.x0, self.y0, self.t0, self.x1, self.y1, self.t1, + self.norad, self.catalog, self.state) class fourframe: @@ -95,15 +95,15 @@ class fourframe: self.zavg, self.zstd, self.zmax, self.znum = hdu[0].data # Generate sigma frame - self.zsig = (self.zmax-self.zavg)/(self.zstd+1e-9) + self.zsig = (self.zmax - self.zavg) / (self.zstd + 1e-9) # Frame properties self.ny, self.nx = self.zavg.shape self.nz = hdu[0].header['NFRAMES'] # Read frame time oselfsets - self.dt = np.array([hdu[0].header['DT%04d' % i] - for i in range(self.nz)]) + self.dt = np.array( + [hdu[0].header['DT%04d' % i] for i in range(self.nz)]) # Read header self.mjd = hdu[0].header['MJD-OBS'] @@ -114,30 +114,29 @@ class fourframe: self.fname = fname # Astrometry keywords - self.crpix = np.array([hdu[0].header['CRPIX1'], - hdu[0].header['CRPIX2']]) - self.crval = np.array([hdu[0].header['CRVAL1'], - hdu[0].header['CRVAL2']]) - self.cd = np.array([[hdu[0].header['CD1_1'], - hdu[0].header['CD1_2']], - [hdu[0].header['CD2_1'], - hdu[0].header['CD2_2']]]) + self.crpix = np.array( + [hdu[0].header['CRPIX1'], hdu[0].header['CRPIX2']]) + self.crval = np.array( + [hdu[0].header['CRVAL1'], hdu[0].header['CRVAL2']]) + self.cd = np.array( + [[hdu[0].header['CD1_1'], hdu[0].header['CD1_2']], + [hdu[0].header['CD2_1'], hdu[0].header['CD2_2']]]) self.ctype = [hdu[0].header['CTYPE1'], hdu[0].header['CTYPE2']] self.cunit = [hdu[0].header['CUNIT1'], hdu[0].header['CUNIT2']] - self.crres = np.array([hdu[0].header['CRRES1'], - hdu[0].header['CRRES2']]) + self.crres = np.array( + [hdu[0].header['CRRES1'], hdu[0].header['CRRES2']]) hdu.close() - + # Compute image properties - self.sx = np.sqrt(self.cd[0, 0]**2+self.cd[1, 0]**2) - self.sy = np.sqrt(self.cd[0, 1]**2+self.cd[1, 1]**2) - self.wx = self.nx*self.sx - self.wy = self.ny*self.sy - self.zmaxmin = np.mean(self.zmax)-2.0*np.std(self.zmax) - self.zmaxmax = np.mean(self.zmax)+6.0*np.std(self.zmax) - self.zavgmin = np.mean(self.zavg)-2.0*np.std(self.zavg) - self.zavgmax = np.mean(self.zavg)+6.0*np.std(self.zavg) + self.sx = np.sqrt(self.cd[0, 0]**2 + self.cd[1, 0]**2) + self.sy = np.sqrt(self.cd[0, 1]**2 + self.cd[1, 1]**2) + self.wx = self.nx * self.sx + self.wy = self.ny * self.sy + self.zmaxmin = np.mean(self.zmax) - 2.0 * np.std(self.zmax) + self.zmaxmax = np.mean(self.zmax) + 6.0 * np.std(self.zmax) + self.zavgmin = np.mean(self.zavg) - 2.0 * np.std(self.zavg) + self.zavgmax = np.mean(self.zavg) + 6.0 * np.std(self.zavg) # Setup WCS self.w = wcs.WCS(naxis=2) @@ -163,10 +162,10 @@ class fourframe: def selection_mask(self, sigma, zstd): """Create a selection mask""" c1 = ndimage.uniform_filter(self.znum, 3, mode='constant') - c2 = ndimage.uniform_filter(self.znum*self.znum, 3, mode='constant') + c2 = ndimage.uniform_filter(self.znum * self.znum, 3, mode='constant') # Add epsilon to keep square root positive - z = np.sqrt(c2-c1*c1+1e-9) + z = np.sqrt(c2 - c1 * c1 + 1e-9) # Standard deviation mask c = z < zstd @@ -177,7 +176,7 @@ class fourframe: c = self.zsig < sigma m2 = np.zeros_like(self.zavg) m2[~c] = 1.0 - self.zsel = m1*m2 + self.zsel = m1 * m2 # Generate points c = self.zsel == 1.0 @@ -189,12 +188,17 @@ class fourframe: return x, y, inum, t, sig - def significant_pixels_along_track(self, sigma, x0, y0, - dxdt, dydt, rmin=10.0): + def significant_pixels_along_track(self, + sigma, + x0, + y0, + dxdt, + dydt, + rmin=10.0): """Extract significant pixels along a track""" # Generate sigma frame - zsig = (self.zmax-self.zavg)/(self.zstd+1e-9) + zsig = (self.zmax - self.zavg) / (self.zstd + 1e-9) # Select c = (zsig > sigma) @@ -207,9 +211,9 @@ class fourframe: t = np.array([self.dt[i] for i in inum]) # Predicted positions - xr = x0+dxdt*t - yr = y0+dydt*t - r = np.sqrt((x-xr)**2+(y-yr)**2) + xr = x0 + dxdt * t + yr = y0 + dydt * t + r = np.sqrt((x - xr)**2 + (y - yr)**2) c = r < rmin return x[c], y[c], t[c], sig[c] @@ -218,7 +222,7 @@ class fourframe: """Extract significant pixels""" # Generate sigma frame - zsig = (self.zmax-self.zavg)/(self.zstd+1e-9) + zsig = (self.zmax - self.zavg) / (self.zstd + 1e-9) # Select c = (zsig > sigma) @@ -239,8 +243,8 @@ class fourframe: # Loop over frames for i in range(self.nz): - dx = int(np.round(dxdt*(self.dt[i]-tref))) - dy = int(np.round(dydt*(self.dt[i]-tref))) + dx = int(np.round(dxdt * (self.dt[i] - tref))) + dy = int(np.round(dydt * (self.dt[i] - tref))) # Skip if shift larger than image if np.abs(dx) >= self.nx: @@ -250,28 +254,23 @@ class fourframe: # Extract range if dx >= 0: - i1min, i1max = dx, self.nx-1 - i2min, i2max = 0, self.nx-dx-1 + i1min, i1max = dx, self.nx - 1 + i2min, i2max = 0, self.nx - dx - 1 else: - i1min, i1max = 0, self.nx+dx-1 - i2min, i2max = -dx, self.nx-1 + i1min, i1max = 0, self.nx + dx - 1 + i2min, i2max = -dx, self.nx - 1 if dy >= 0: - j1min, j1max = dy, self.ny-1 - j2min, j2max = 0, self.ny-dy-1 + j1min, j1max = dy, self.ny - 1 + j2min, j2max = 0, self.ny - dy - 1 else: - j1min, j1max = 0, self.ny+dy-1 - j2min, j2max = -dy, self.ny-1 + j1min, j1max = 0, self.ny + dy - 1 + j2min, j2max = -dy, self.ny - 1 zsel = np.where(self.znum == i, self.zmax, 0.0) ztrk[j2min:j2max, i2min:i2max] += zsel[j1min:j1max, i1min:i1max] return ztrk def __repr__(self): - return "%s %dx%dx%d %s %.3f %d %s" % (self.fname, - self.nx, - self.ny, - self.nz, - self.nfd, - self.texp, - self.site_id, - self.observer) + return "%s %dx%dx%d %s %.3f %d %s" % (self.fname, self.nx, self.ny, + self.nz, self.nfd, self.texp, + self.site_id, self.observer) diff --git a/stvid/utils.py b/stvid/utils.py index 865c20d..2ef6ae7 100644 --- a/stvid/utils.py +++ b/stvid/utils.py @@ -13,8 +13,8 @@ def get_sunset_and_sunrise(tnow, loc, refalt_set, refalt_rise): # Get time nmjd = 64 mjd0 = np.floor(tnow.mjd) - mnow = tnow.mjd-mjd0 - mjd = np.linspace(mjd0-1.0, mjd0+3.0, nmjd) + mnow = tnow.mjd - mjd0 + mjd = np.linspace(mjd0 - 1.0, mjd0 + 3.0, nmjd) t = Time(mjd, format='mjd', scale='utc') # Get sun position @@ -28,8 +28,10 @@ def get_sunset_and_sunrise(tnow, loc, refalt_set, refalt_rise): # Compute altitude extrema de = np.mean(pos.dec) - minalt = np.arcsin(np.sin(loc.lat)*np.sin(de)-np.cos(loc.lat)*np.cos(de)) - maxalt = np.arcsin(np.sin(loc.lat)*np.sin(de)+np.cos(loc.lat)*np.cos(de)) + minalt = np.arcsin( + np.sin(loc.lat) * np.sin(de) - np.cos(loc.lat) * np.cos(de)) + maxalt = np.arcsin( + np.sin(loc.lat) * np.sin(de) + np.cos(loc.lat) * np.cos(de)) # Never sets, never rises? if minalt > min(refalt_set, refalt_rise): @@ -38,55 +40,48 @@ def get_sunset_and_sunrise(tnow, loc, refalt_set, refalt_rise): return "sun never rises", t[0], t[0] # Prevent discontinuities in right ascension - dra = pos[-1].ra-pos[0].ra + dra = pos[-1].ra - pos[0].ra ra = pos.ra.degree if dra < 0.0: c = pos.ra.degree < 180.0 - ra[c] = pos[c].ra.degree+360.0 + ra[c] = pos[c].ra.degree + 360.0 # Set up interpolating function fra = interpolate.interp1d(t.mjd, ra) fde = interpolate.interp1d(t.mjd, pos.dec.degree) # Get GMST - gmst0 = Time(mjd0, - format='mjd', + gmst0 = Time(mjd0, format='mjd', scale='utc').sidereal_time('mean', 'greenwich') # Get transit time - mtransit = np.mod((fra(mjd0)*u.degree-loc.lon-gmst0)/(360.0*u.degree), 1.0) + mtransit = np.mod( + (fra(mjd0) * u.degree - loc.lon - gmst0) / (360.0 * u.degree), 1.0) while True: - gmst = gmst0+360.985647*u.deg*mtransit - ra = fra(mjd0+mtransit)*u.deg - ha = gmst+loc.lon-ra - mtransit -= ha/(360.0*u.deg) + gmst = gmst0 + 360.985647 * u.deg * mtransit + ra = fra(mjd0 + mtransit) * u.deg + ha = gmst + loc.lon - ra + mtransit -= ha / (360.0 * u.deg) if np.abs(ha.degree) < 1e-9: break # Hour angle offset - ha0 = np.arccos((np.sin(refalt_set) - - np.sin(loc.lat) - * np.sin(np.mean(pos.dec))) - / (np.cos(loc.lat) - * np.cos(np.mean(pos.dec)))) - + ha0 = np.arccos( + (np.sin(refalt_set) - np.sin(loc.lat) * np.sin(np.mean(pos.dec))) / + (np.cos(loc.lat) * np.cos(np.mean(pos.dec)))) + # Get set time - mset = mtransit+ha0/(360.0*u.deg) + mset = mtransit + ha0 / (360.0 * u.deg) while True: - gmst = gmst0+360.985647*u.deg*mset - ra = fra(mjd0+mset)*u.deg - de = fde(mjd0+mset)*u.deg - ha = gmst+loc.lon-ra - alt = np.arcsin(np.sin(loc.lat) - * np.sin(de) - + np.cos(loc.lat) - * np.cos(de) - * np.cos(ha)) - dm = (alt-refalt_set)/(360.0 - * u.deg - * np.cos(de) - * np.cos(loc.lat) - * np.sin(ha)) + gmst = gmst0 + 360.985647 * u.deg * mset + ra = fra(mjd0 + mset) * u.deg + de = fde(mjd0 + mset) * u.deg + ha = gmst + loc.lon - ra + alt = np.arcsin( + np.sin(loc.lat) * np.sin(de) + + np.cos(loc.lat) * np.cos(de) * np.cos(ha)) + dm = (alt - refalt_set) / (360.0 * u.deg * np.cos(de) * + np.cos(loc.lat) * np.sin(ha)) mset += dm # Break loop or find sunset on next day @@ -97,25 +92,20 @@ def get_sunset_and_sunrise(tnow, loc, refalt_set, refalt_rise): mset += 1.0 # Set set time - tset = Time(mjd0+mset.value, format='mjd', scale='utc') + tset = Time(mjd0 + mset.value, format='mjd', scale='utc') # Get rise time - mrise = mtransit-ha0/(360.0*u.deg) + mrise = mtransit - ha0 / (360.0 * u.deg) while True: - gmst = gmst0+360.985647*u.deg*mrise - ra = fra(mjd0+mrise)*u.deg - de = fde(mjd0+mrise)*u.deg - ha = gmst+loc.lon-ra - alt = np.arcsin(np.sin(loc.lat) - * np.sin(de) - + np.cos(loc.lat) - * np.cos(de) - * np.cos(ha)) - dm = (alt-refalt_rise)/(360.0 - * u.deg - * np.cos(de) - * np.cos(loc.lat) - * np.sin(ha)) + gmst = gmst0 + 360.985647 * u.deg * mrise + ra = fra(mjd0 + mrise) * u.deg + de = fde(mjd0 + mrise) * u.deg + ha = gmst + loc.lon - ra + alt = np.arcsin( + np.sin(loc.lat) * np.sin(de) + + np.cos(loc.lat) * np.cos(de) * np.cos(ha)) + dm = (alt - refalt_rise) / (360.0 * u.deg * np.cos(de) * + np.cos(loc.lat) * np.sin(ha)) mrise += dm # Break loop or find sunrise on next day @@ -126,6 +116,6 @@ def get_sunset_and_sunrise(tnow, loc, refalt_set, refalt_rise): mrise += 1.0 # Set rise time - trise = Time(mjd0+mrise.value, format='mjd', scale='utc') + trise = Time(mjd0 + mrise.value, format='mjd', scale='utc') return "sun rises and sets", tset, trise diff --git a/update_tle.py b/update_tle.py index 9909b17..63d48eb 100755 --- a/update_tle.py +++ b/update_tle.py @@ -3,7 +3,6 @@ from __future__ import print_function import configparser import argparse from spacetrack import SpaceTrackClient -import subprocess from shutil import copyfile import datetime from io import BytesIO