import os import sys import spiceypy as spice """ Usage: python3 nh_spin_state.py \ [--outfile=nh_spin.asc] \ [--csvfile=nh_spin.csv] \ tsoc/recon/ck/merged_nhpc_20*.bc \ tsoc/common/lsk/naif0010.tls \ tsoc/common/sclk/new-horizons_0664.tsc \ tsoc/common/fk/nh_v220.tf """ # BLE: PROBABLY RE-ADD THESE TO FIX OUTPUT? #if not ("SPINDEBUG" in os.environ): # spice.errprt("set","none") # spice.erract("set","report") # spice.errdev("set","xyz") ### Parse input arguments ### ### --outfile=nh_spin.asc name of file to receive the start or end time ### of each contiguous period of time at one state (spinning or 3-axis) ### written on separate lines ### ### --csvfile=nh_spin.csv name of file to receive the start and end ### time of each contiguous period of time at one state (spinning or ### 3-axis) written on one line ### ### All other arguments are assumed to be input SPICE Kernels ### - Need C-Kernels, an SCL-Kernel and an LS-Kernel ### ### WARNING: The time range is hard-coded below. This line should be ### changed to extend the output range: ### while utc<'2023': outfile = False ### Initailze outputs to go to STDOUT csvfile = False ### " ### Loop over command-line arguments for fn in sys.argv[1:]: ### file for one start or end time per line if fn[:10]=='--outfile=': outfile=fn[10:] ### file for one start and end time per line elif fn[:10]=='--csvfile=': csvfile=fn[10:] else: ### argument is a SPICE kernel try: spice.furnsh(fn) except: pass ### Start before launch utc='2006-01-01T12:00:00' ### Convert UTC to ET, set linear-search timestep et=spice.utc2et(utc) deltad = 1024 ### Conversion from radian/s to RPM ### radian/s x DPR degree/radian x 60 s/Minute x 1/360 Revolution/degree = RPM k = spice.dpr() * 60. / 360. ######################################################################## ### Method to set attitude mode and calculate rotation rate around Y ### - Gets state matrix from C-kernels at time ET (SXFORM) ### - Convert state matrix to angular velocities around Y (XF2ARAV) ### - Use 2.5RPM as knife-edge trigger to set attitude mode ### - Exception is assumed to be a SPICE call failure; set mode to missing def getrpm(et): spice.reset() try: state = spice.sxform('NH_SPACECRAFT','J2000', et) rot,av = spice.xf2rav(state) rpm = spice.vnorm(av)*k if rpm > 2.5: s='Spinning' else : s='3-axis' return rpm,s,False except: if "SPINDEBUG" in os.environ: import traceback traceback.print_exc() print( utc ) return -1e32,'Missing',True ######################################################################## def binsearch(et0Arg,s0,et1Arg): """ Method to do do a binary search in time to find a mode transition to the nearest second binsearch(et0, s0, et1) et0 - ephemeris time at initial state s0 - state at et0 et1 - ephemeris time at second, different state - Bidirectional; et0 may be less than or greater than et1 - Assumes abs(et0-et1) is an integral power of two - Assumes the state at et1 is not s0 - Finds a time less than 1s closer to et0 than a mode transition away from s0 - Treats the two possible modes other than s0 as "not s0" - Weaknesses - Possible resolution and/or overflow problem; unlikely with double precision ETs ca. 21st century - Due to the possibility of gaps in coverage and dependent on the initial time separation between et0 and et1, this is not rigorous, but if used from either end, the transition times will not overlap """ ### Copy endpoint times to local variables et0,et1 = et0Arg, et1Arg ### Do binary search to resolution of 1s while abs(et1-et0)>1.1: ### Possible infinite loop etMid = (et0+et1) / 2.0 ### Midpoint time ignore,s,m = getrpm(etMid) ### Get state at midpoint if s==s0: et0 = etMid ### If state is not s0 set et0 to midpoint else: et1 = etMid ### Else set et1 to midpoint sys.stdout.flush() return et0 ### Return last time on et0 side of transition ### String to indicate if binary search result is at a gap in coverage mdf = lambda missing,fp: ( ' (At least 1 coverage gap%s)' % (fp,)) if missing else '' ### Initialize variables s (state) and missing to missing and True ignore,s,missing = getrpm(None) lastet = et prevLine = False ### Stack to keep track of forward-in-time transitions to or from ### Spinning or 3-axis mode stack = [] stacklasts = '' ######################################################################## ### 2014: currently set for NH SPICE Archive Release 0006; C-Kernel ### coverage ends in December, 2020. while utc<'2023': ### Save last s and missing values, get new values lasts,lastmissing = s,missing rpm,s,missing = getrpm(et) ### Check whether to do binary search for transitions: if s!=lasts: if not lastmissing: ### From lastet to et if state at lastet is not missing et0 = binsearch(lastet,lasts,et) ### - rpm0,s0,missing0 = getrpm(et0) utc0 = spice.et2utc(et0,'ISOC',0) rpm1,s1,missing1 = getrpm(et0+1) ### Put state, utc, rpm, and whether transition is to missing ### state onto stack as end time of this state ### Keep track, in n, of the number of contiguous times for this ### transition (e.g. Spinning, missing, Spinning, missing, ...) if s0==stacklasts: n = n+1 else : n,stacklasts = 0, s0 stack.append( [s0, '%s %8.3f %s end%s' % (utc0,rpm0,s0, mdf(missing1,' following'),),n] ) if not missing: ### From et to lastet if state at et is not missing et0 = binsearch(et,s,lastet) rpm0,s0,missing0 = getrpm(et0) utc0 = spice.et2utc(et0,'ISOC',0) rpm1,s1,missing1 = getrpm(et0-1) ### Put state, utc, rpm, and whether transition is from missing ### state onto stack as end time of this state ### Keep track, in n, of the number of contiguous times for this ### transition (e.g. Spinning, missing, Spinning, missing, ...) if s0==stacklasts: n = n+1 else : n,stacklasts = 0, s0 stack.append( [s0, '%s %8.3f %s start%s' % (utc0,rpm0,s0, mdf(missing1,' preceding'),),n] ) ###sys.stdout.flush() ### Save et to lastet, step et, determine utc for while utc<2014: above lastet = et et += deltad utc = spice.et2utc(et,'ISOC',0) ######################################################################## ### Work backwards down stack built above, collapsing start and ### end of contiguous sequence of same-state transitions with ### missing state each into one item on a new stack, stack2 stack2 = [] lasts = '' while stack: ### Pop state, line, and n (count) off of stack ### - n will be odd when a new state (s !== lasts) is popped ### - n will be 0 when the last of a contiguous sequence of the same ### state is popped s,line,n = stack.pop() ### Skip subsequent contiguous transitions with the same state ### and a non-zero count; the zero count is the initial transition if s==lasts and n>0: continue assert (n&1) or n==0 lasts = s ### Determine intra-sequence gap comments based on n ### - n==0 or 1, no comment if n<2 : sfx='' ### - n==3: one gap was found elif n==3: sfx=' (At least 1 coverage gap within)' ### - n==5,7,...: more than one gap was found else : sfx=' (At least %d coverage gaps within)' % (n>>1,) ### Append to stack2, which will have pairs of alternating start/end ### elements, with each successive pair alternating Spinning and 3-axis ### N.B. stack2 is reversed in time (last item is first time) stack2.append( '%s%s' % (line,sfx,) ) ######################################################################## ### Open --outfile-... file for stack2 output of one transition per line ### BLE: Changed 'wb' to 'w' to migrate to python3.5 . if outfile: f = open(outfile,'w') else : f = sys.stdout ### Copy stack2 to temporary stack stack2tmp = stack2[:] ### Print items from stack2, then close file for s2 in stack2: f.write( '%s\n' % (stack2tmp.pop(),) ) if outfile: f.close() ######################################################################## ### Open --csvfile-... file for stack2 output of two transitions per line ### BLE: Changed 'wb' to 'w' to migrate to python3.5 . if csvfile: f = open(csvfile,'w') else : f = sys.stdout ### Write header f.write( 'Start UTC,Stop Utc,RPM at entry,RPM at exit,State,' + 'Gaps comment 1,Gaps comment 2,Gaps comment 3\n' ) ### Lambda functions to reformat pairs of lines ### - pop item from stack, remove closing parens, leaving opening parens get1 = lambda stck: stck.pop().strip(' )\n' ).replace(') (','(' ).replace(')(','(' ) ### - split at opening parens, return items after first paren extras = lambda line: line.split('(')[1:] while stack2: ### Get next pair of [start,end] transitions try: line1,line2 = get1(stack2), get1(stack2) except IndexError: break ### Extract space-separated times, rpms, states utc1,rpm1,state1 = line1.split()[:3] utc2,rpm2,state2 = line2.split()[:3] ### Ensure states within pair are the same assert state1==state2 ### Reformat any comments, ensure there are at least 3, keep first 3 gaps = (extras(line1) + extras(line2) + ['','','',''])[:3] ### Write line as Comma-Separated Values (CSVs) f.write( '%s\n' % ','.join( [utc1,utc2,state1[:4],rpm1,rpm2]+gaps ) ) ### Close file if csvfile: f.close()