#!/usr/local/bin/python3.11 #simple/limited implementation of non-sidereal tracking # will move fiducial on active guider to compensate # for non-sidereal motion. LIMITED!! # 27-Dec 2023 SPJ, started working from pseudo-code import numpy as np import scipy #from astropy.io import fits from astropy.io import ascii #import matplotlib #matplotlib.use('Agg') #import matplotlib.pyplot as plt #from matplotlib.ticker import MultipleLocator, FormatStrFormatter #import itertools #import random import datetime import time import argparse #import importlib import os, sys from os.path import exists #import re #import glob from tcs_lib import tcs_event from hetlib import tcsutils as tu from TCSNamed import TCSNamed import tcssubsystem as ts def read_eph(eph_file): #read in the ephemerides file from JPL HORIZONS print('reading in file: '+str(eph_file[0])) indat = ascii.read(eph_file[0]) #print(indat) #convert dates into datetime objects: dt = np.array([datetime.datetime.strptime(dd+' '+tt, '%Y-%b-%d %H:%M') for dd,tt in zip(indat['col1'],indat['col2'])]) #print(dt[0]) #convert RA and Dec into decimal degrees I guess ra = np.array([15*(hh+mm/60+ss/3600) for hh,mm,ss in zip(indat['col3'],indat['col4'],indat['col5'])]) dec = np.array([dd+mm/60+ss/3600 for dd,mm,ss in zip(indat['col6'],indat['col7'],indat['col8'])]) #print(ra[0],dec[0]) return(dt,ra,dec) def parse_input(): myDescription = """ Offset guider fiducial to accomplish non-sidereal tracking (over limited range!) Ephemeris file must have format like this on each line, with UT date/time, RA, and Dec (ICRF): 2023-Nov-29 00:00:00 02 20 15.15 +12 33 52.3 Examples: ./nst.py amalthea_20231129.txt (use ephemerides file listed and default 1min updates) ./nst.py amalthea_20231129.txt --step 1.5 (use 1.5 minute time step between fiducal offsets instead of default 1 min) """ parser = argparse.ArgumentParser(description=myDescription,formatter_class=argparse.RawDescriptionHelpFormatter) parser.add_argument('eph_file',metavar='eph_file', type=str, nargs=1, help='ASCII file with ephemerides from HORIZONS: Date, RA, DEC') parser.add_argument('--step', metavar='step', type=float, nargs='?', default=1, help='time step (minutes) between fiducial offsets') global args #cheating args = parser.parse_args() #sanity check... not a crazy number? if args.step > 30.0: print('you have requested updates every '+str(args.step)+' minutes, which is unreasonably slow') print(' please try again with a different step size') sys.exit() if args.step < 0.2: print('you have requested updates every '+str(args.step)+' minutes, which is unreasonably fast') print(' please try again with a different step size') sys.exit() return(args) def check_continue(): #check if anyone has created a "stop" file if (exists('stop')): print(' You have asked that this script be stopped! it will stop now...') os.remove('stop') sys.exit() return False else: return True def do_nst(dt,ra,dec,gact,g,step): print('') print(' !!!! non-sidereal tracking activated !!!! ') print('') print(' will make fiducial corrections every '+str(step)+' minute(s)') print(' (touch a file called stop in this directory to stop this script) ') print('') check_continue() print('activated at: '+str(gact)) next_update = gact+datetime.timedelta(minutes=step) print('next update: '+str(next_update)) #interpolate to get exact position at activation time # first, convert timestamps to seconds for easier interpolation dts = np.array([(t - dt[0]).total_seconds() for t in dt]) #print(dts) #print(ra) #then convert guider activation time into seconds as well gacts = (gact-dt[0]).total_seconds() #now do the interpolation: #ra0,dec0 = np.interp(gacts, dts, ra), np.interp(gacts, dts, dec) #higher order interpolation: ra0 = scipy.interpolate.interp1d(dts, ra, kind='cubic')(gacts) dec0 = scipy.interpolate.interp1d(dts, dec, kind='cubic')(gacts) #print(ra0,dec0) #print(_ra0,_dec0) #set the "last" position to this one ra_last,dec_last = ra0,dec0 keeplooping=True while(keeplooping): #infinite loop check_continue() #stop if user says stop #is it time for the next update? if datetime.datetime.now()>=next_update: #interpolate to get exact position at next_update time n = datetime.datetime.now() dtsnow = (n-dt[0]).total_seconds() #ra_now,dec_now = np.interp(dtsnow, dts, ra), np.interp(dtsnow, dts, dec) #higher order interpolation: ra_now = scipy.interpolate.interp1d(dts, ra, kind='quadratic')(dtsnow) dec_now = scipy.interpolate.interp1d(dts, dec, kind='quadratic')(dtsnow) #print(ra_now,dec_now) #compare with last position (VERIFY that this sign is correct??) dra,ddec = 3600*(ra_now - ra_last), 3600*(dec_now - dec_last) #print(dra,ddec) print(' updating g'+str(g)+' fiducial by dra,ddec='+str(dra)+','+str(ddec)) #send this as fiducial offset with tcs: if (g==1): #print('g1 offset') tcs.Guider1_offset_fiducial(dra_asec=dra,ddec_asec=ddec,compensate="true") #tcs.() else: #print('g2 offset') tcs.Guider2_offset_fiducial(dra_asec=dra,ddec_asec=ddec,compensate="true") #tcs.() ra_last,dec_last=ra_now,dec_now next_update+=datetime.timedelta(minutes=step) print('next update: '+str(next_update)) time.sleep(0.1) #wait this many seconds before checking again #print(tcs) return None if __name__ == "__main__": if(exists('stop')): os.remove('stop') #so it won't immediately stop! #parse input: args=parse_input() #print(args) #get ephemerides dt,ra,dec = read_eph(args.eph_file) #print(dt,ra,dec) #ra/dec are in decimal degrees) #do sanity checks on those - current time OK? enough time remaining? # current time? n = datetime.datetime.now() #print(n) #check if current time is within the range of times mi,ma = dt.min(),dt.max() #print(mi,ma) if n>mi and n