PRO PMTRAS_DBASE_GENERATOR, initial_time, total_duration, $ DBSOLUTION=dbsolution, SAVEFILE=savefile, $ NOLOG=nolog, PMTRAS_DIAGNOSTIC=diag_plots, $ SHORT_INTERVALS=short_intervals, $ SC_SUN_OFFSET= sc_sun_offset ; ; Generates potential elements of LOCAL_ROLL_DBASE.dat ; Data points will be at all values of s/c clock that are within obs_time_interval AND that ; are integral multiples of 64s. ; Multiple calls to pmtras_analysis (about one per orbit) are used. ; ; ; INPUTS: initial_time = start time of analysis interval in anytime format. eg ('yyyy-mmm-dd hh:mm') ; total_duration = total duration of analysis interval (minutes) ; ; KEYWORDS: ; DBSOLUTION = output structure containing the solution. ; SAVEFILE = Name of .sav file in which to save solution structure, dbsolution. ; Typical file name is of form 'PMTRASyyyymmmdd-dd.SAV' ; Default file name is 'PMTRAS_DBASE_GEN.SAV' ; NOLOG = Set to suppress output to 'PMTRAS_GEN_LOG.txt'. Default = 0. ; PMTRAS_DIAGNOSTIC = pmtras_diagnostic code. (Default=128.) ; SHORT_INTERVALS = Set to use 32-minute intervals with 16-minute offset for pmtras analysis intervals. ; Default is 96-minute intervals with 0-minute offset. ; SC_SUN_OFFSET = 2-element vector of [+W, +N] sc offset with respect to Sun (degrees). Needed if ; offset is >~ 1 degree. Default = [0,0]. ; ; USAGE EXAMPLE: PMTRAS_DBASE_GENERATOR, '2002-apr-21 00:30', 1440 ; generates 24-hour solution. ; ; CALLS: anytim ; hsi_any2sctime ; hsi_packet ; hsi_pmtras_multifit ; hsi_sctime2any ; pmtras_analysis ; ; 9-28-02 gh Initial version - local database parameters incomplete. ghurford@ssl.berkeley.edu ; 10-01-02 gh Local database parameters complete. ; 11-08-02 gh Minor correction to time_extension and reltimeout ; 11-12-02 gh Clean up assignment of roll_quality ; Add provision for log generation. ; 11-13-02 gh Modify input to support multi-orbit analyses ; Add provision for SAVE file generation ; Add keyword for pmtras_diagnostic level, with default being a final phase-residual plot. ; 12-31-02 gh Update call to pmtras_analysis. ; 13-jan-03 gh Fix bug which caused crash if 1st orbit had no solution. ; 5-May-03 gh Added line count message to screen output. ; 16-Jul-03 gh Adapt to modifications in hsi_pmtras_linfit calling arguments. ; Exclude outliers from list of stars in dbsoln. ; 22-Jul-03 gh Add PMTRAS_DIAGNOSTIC keyword. ; Correct star parameter bug introduced 16-July ; Debug indexing errors in star-specific tags. ; 29-Jul-03 gh Add provision for case where there is too little data in one orbit to process. ; gh Set acode==2 ; 30-Jul-03 gh Test that pmtras_analysis returned a structure for this orbit. ; 1-Aug-03 gh Correct bug in 29-July change ; 4-Aug-03 gh Add /NOSTOP keyword to pmtras_analysis call to force continuation after a failed analysis. ; 12-Sep-03 gh Use hsi_pmtras_multifit to support option of quadratic fits within 64s intervals ; Require est_error > 0 for good (AA=11) data quality ; Set acode==4 ; 15-Oct-03 gh Add SHORT_INTERVALS keyword to support verification of star identifiation. ; 9-Jun-04 gh Added SC_SUN_OFFSET keyword. ; 2-Jun-06 gh Open a new 'pmtras_gen_log.txt' file if it does not exist. ; 8-Jun-07 gh Debug last year's 'improvement'. ; ; Define 1-byte analysis code which flags program version ; LSbit will be used to flag 'manual edits'. ; ; acode = 0B ; until 28-July-03 Initial version ; acode = 2B ; 31-July-03 to 11-Sep-03 Outliers excluded from dbsoln.stars; extra_time_offset=0 acode = 4B ; 12-Sep-03 Optional quad. interpolation; require est_error>0 for good quality ; ; ; Define some fixed parameters ; TWOPI = 2.D * !DPI TWOPIMR = TWOPI * 1.D6 ; microradians per rotation IF KEYWORD_SET(short_intervals) THEN BEGIN ; data points will correspond to integer*orbit_time + orbit_offset orbit_time = 32. orbit_offset = 16. ENDIF ELSE BEGIN orbit_time = 96. ; nominal interval for one pmtras_analysis (minutes) orbit_offset = 0. ; ENDELSE time_extension = [-96, +32] ; required to ensure data for +-32s of all output points. IF N_ELEMENTS(savefile) EQ 0 THEN savefile = 'PMTRAS_DBASE_GEN.SAV' IF N_ELEMENTS(diag_plots) EQ 0 THEN diag_plots = 128 IF N_ELEMENTS(sc_sun_offset) EQ 0 THEN sc_sun_offset = [0,0] ; ; Open log file and define packet object. ; fileinfo = FILE_INFO('PMTRAS_GEN_LOG.txt') IF fileinfo.exists NE 0 THEN OPENU, dbglogfile, 'PMTRAS_GEN_LOG.txt', /GET_LUN, /APPEND $ ELSE OPENW, dbglogfile, 'PMTRAS_GEN_LOG.txt', /GET_LUN nlines = 0 opkt = hsi_packet() ; ; Use PMTRAS_ANALYSIS to generate a solution in structure, dbsoln ; norb = ROUND(total_duration / orbit_time) > 1 ; number of pmtras_analysis runs to be done duration = FLOAT(total_duration) / norb ; duration of each pmtras_analysis analysis run ; ; Begin loop over orbits. ; FOR nor = 0, norb-1 DO BEGIN toff = 60 * (duration * [nor, nor+1] + orbit_offset) ; vector defining time(minutes) within day for 1 analysis obs_time_interval = anytim(initial_time) + toff opkt -> set, app_id=154, obs_time_interval=obs_time_interval+ time_extension pmtras_data = opkt->getdata() IF SIZE(pmtras_data, /TYPE) EQ 2 THEN CONTINUE ; no data in this orbit - CONTINUE at next orbit pmtras_analysis, packet_object=opkt,pmtras_diagnostic=diag_plots, dbase_solution=dbsoln, $ SC_SUN_OFFSET=sc_sun_offset, /NOLOOKUP, /NOSTOP IF SIZE(dbsoln, /TYPE) EQ 2 THEN CONTINUE ; no solutions in this orbit - CONTINUE at next orbit ; ; Define a series of s/c reference times which are integral multiples of 64s. ; These output times are expressed in array, reltimeout, as times relative to t0_ref. ; nblip = N_ELEMENTS(dbsoln.reltime) ; number of blips in dbsoln scobstime = hsi_any2sctime(obs_time_interval) ; obs_time_interval as a s/c clock structure obs_time = DOUBLE(scobstime.seconds) + FLOAT(scobstime.bmicro)/2L^20 ; obs_time_interval as s/c clock (seconds) tinref = DOUBLE(dbsoln.t0_ref.seconds) + FLOAT(dbsoln.t0_ref.bmicro)/2L^20 ; time is s/c clock (seconds) tinfirst = tinref + dbsoln.reltime[0] ; time of 1st solution datum tinlast = tinref + dbsoln.reltime[nblip-1] ; time of last solution datum toutfirst = (ROUND(tinfirst/64) +1) * 64 ; 1st scref time -- has input data for +- 32s toutlast = (ROUND(tinlast /64) -1) * 64 ; last scref time -- has input data for +- 32s ndt0 = FIX((obs_time[0] - toutfirst)/64 +1) > 0 ndt1 = FIX((toutlast - obs_time[1])/64 +1) > 0 toutfirst = toutfirst + ndt0*64 ; forces toutfirst to be after obs_time_interval[0] toutlast = toutlast - ndt1*64 ; forces toutlast to be before obs_time_interval[1] ntimes = ROUND((toutlast - toutfirst)/64) + 1 ; number of output points IF ntimes LE 0 THEN CONTINUE ; not enough data - CONTINUE at next orbit reltimeout = toutfirst - tinref + INDGEN(ntimes)*64 ; array of scref times as 'reltimes' ; ; Define the 1-orbit output structure, dbsol, and populate the sctime values. ; star = {detections: 0, $ ; detections of this star in this interval id: 0, $ ; Harvard ref# of star intensity: 0, $ ; average blip total av_offset: 0.0, $ ; average offset from overall fit in this interval rms_scatter: 0.0 } ; rms scatter, relative to av_offset (defined for ndet>1) stars = REPLICATE (star,10) ; can support up to 10 stars dbsol = {sctime: 0L, $ ; in s/c clock seconds, should be a multiple of 64 roll_phase: 0L, $ ; in microradians approx_angvel: 0B, $ ; approximate angular velocity, in radians/second * 2^7 roll_quality: 0B, $ ; roll_quality byte blipcount: 0, $ ; number of blips in this interval starcount: 0B, $ ; number of unique stars in this interval roll_period: 0.0, $ ; exact roll period (s) est_error: 0.0, $ ; estimated error in roll_phase (radians) analysis_code: 0B, $ ; TBD code to indicated software configuration/parameters stars: stars } dbsol.analysis_code = acode dbsol = REPLICATE (dbsol, ntimes) dbsol.sctime = toutfirst + INDGEN(ntimes)*64 ; ; Begin loop over output points, processing 64 seconds of data at a time. ; Details of processing depend on whether 0,1,2 or >2 blips were detected. ; FOR n = 0, ntimes-1 DO BEGIN j0 = WHERE (ABS(dbsoln.reltime - reltimeout[n]) LE 32, j0count) ; j0=array of indices to dbsoln in this interval IF j0count GT 1 THEN BEGIN timej0 = dbsoln.reltime[j0] ; integral seconds anglej0 = dbsoln.posn_angle[j0] ; DP monatonically increasing radians fitparm = hsi_pmtras_multifit (timej0, anglej0, est_error, yfitj0, okj) ; linear fit with outlier rejection ; ; Revise j-array, jcount, anglej and yfitj to exclude outliers. ; okj is an array of indices that identify elements of the j0 array of indices that in turn identify the time-selected ; subset of blips. j = j0(okj) ; revise j to use only nonrejected blips jcount = N_ELEMENTS(j) ; j is >= 2. timej = timej0 [okj] anglej = anglej0[okj] yfitj = yfitj0 [okj] ; ; Load tags in dbsol structure timesbracket = timej[0] LT reltimeout[n] AND timej[jcount-1] GT reltimeout[n] ; true if blips phasen = fitparm[0] + fitparm[1]*reltimeout[n] + fitparm[2]*reltimeout[n]^2 ; in radians dbsol[n].roll_phase = LONG((phasen MOD TWOPI) * 1.D6) ; in microradians dbsol[n].roll_period = FLOAT(TWOPI/fitparm[1]) dbsol[n].approx_angvel = BYTE(fitparm[1]*128) dbsol[n].blipcount = jcount dbsol[n].est_error = est_error ; = rms(ok pts) / SQRT(nok-2) ; bracket reltimeout[n] ; ; Determine accuracy bits for roll_quality ; Summary of criteria for accuracy bits, AA: ; 1 arcminute = 0.0003 radians ; AA=11 0 < est_error < 0.0003 AND jcount >2 AND blip times bracket reltimeout ; AA=10 est_error < 0.003 AND jcount >2 AND blip times bracket reltimeout ; AA=01 est_error > 0.003 OR jcount =2 OR blip times do not bracket reltimeout ; AA=00 jcount <2 IF jcount GT 2 AND timesbracket THEN BEGIN IF (est_error LT 0.0003) AND (est_error GT 0) THEN dbsol[n].roll_quality = 216 $ ELSE IF est_error LT 0.003 THEN dbsol[n].roll_quality = 144 $ ELSE dbsol[n].roll_quality = 88 ENDIF ELSE BEGIN IF jcount LT 2 THEN dbsol[n].roll_quality = 8 $ ELSE dbsol[n].roll_quality = 144 ENDELSE ; ; Generate a list of stars detected in the n'th interval hrnj = dbsoln.hrn[j] ampj = dbsoln.intensity[j] paj = dbsoln.posn_angle[j] sorthrnj = hrnj(SORT(hrnj)) hrnlist = sorthrnj(UNIQ (sorthrnj)) ; array of unique HRN values for this interval nstars = N_ELEMENTS(hrnlist) ; number of unique stars for this interval dbsol[n].starcount = nstars ; ; Begin loop over stars detected in the n'th interval ; FOR k=0, nstars-1 DO BEGIN jk = WHERE(hrnj EQ hrnlist[k], jkcount) ; jkcount always >0. diffjk = anglej[jk] - yfitj[jk] dbsol[n].stars[k].detections = jkcount dbsol[n].stars[k].id = hrnlist[k] dbsol[n].stars[k].intensity = MEAN(ampj[jk]) dbsol[n].stars[k].av_offset = MEAN(diffjk) IF jkcount GT 1 THEN BEGIN result = MOMENT(diffjk, SDEV=sdev) dbsol[n].stars[k].rms_scatter = sdev ; std deviation for this star only ENDIF ; rms_scatter defaults to 0 if jkcount=1 ENDFOR ENDIF ELSE BEGIN dbsol[n].roll_quality = 8 ; no PMTRAS blips ENDELSE ENDFOR ; ; Generate one line of log entry ; now = SYSTIME() current_time = STRMID(now,8,2)+'-'+STRMID(now,4,3)+'-'+STRMID(now,20,4)+STRMID(now,10,6) ; 17 char:'yyyy-mmm-dd hh:mm' starttime = hsi_sctime2any(scobstime[0], /STIME) duration = (obs_time[1]-obs_time[0])/60. npt = n_elements(dbsol) dummy = WHERE(dbsol.roll_quality GE 192, ngood) dummy = WHERE(dbsol.roll_quality GE 128, nf) dummy = WHERE(dbsol.roll_quality GE 64, nu) nfair = nf - ngood nuncert = nu - nf nosoln = npt - nu averror = MEAN(dbsol.est_error) maxerror = MAX(dbsol.est_error) IF maxerror GT 0.0003 OR npt NE ngood THEN logflag = '*' ELSE logflag = ' ' IF KEYWORD_SET(nolog) EQ 0 THEN $ PRINTF, dbglogfile, starttime, duration, npt, ngood, nfair, nuncert, nosoln, averror, maxerror, $ current_time, logflag, FORMAT='(A17, F6.1, 5I3, 2F8.5, 1X, A17, A2)' nlines = nlines + 1 IF nor EQ 0 OR N_ELEMENTS(dbsolution) EQ 0 THEN dbsolution = dbsol $ ELSE dbsolution = [dbsolution,dbsol] ; dbsolution holds multiorbit solution ENDFOR SAVE, dbsolution, FILENAME=savefile FREE_LUN, dbglogfile ; ; Generate some summary information ; npts = N_ELEMENTS(dbsolution) PRINT, 'PMTRAS_DBASE_GENERATOR wrote ', npts, ' points to ', savefile IF nlines GT 0 THEN PRINT, 'PMTRAS_DBASE_GENERATOR added ', nlines, ' lines to PMTRAS_GEN_LOG.txt' PLOT, dbsolution.est_error, YTITLE='radians', TITLE='EST_ERRORS' END