PRO PMTRAS_DBASE_PLOT, dbsolution, plotmode
;
; Generates selected plots of PMTRAS database in structure, dbsolution
;
; Set bits in plotmode value select individual plots:
;	Bit 0 = 1 ==> AV_OFFSET vs time						(eg PLOTMODE = 1)
;	Bit 1 = 1 ==> AV_OFFSET histogram					(eg PLOTMODE = 2)
;	Bit 2 = 1 ==> Orbit-to-orbit phase shift vs time 	(eg PLOTMODE = 4)
;	Bit 3 = 1 ==> Roll_period vs time					(eg PLOTMODE = 8)
;	Bit 4 = 1 ==> Est_error vs time						(eg PLOTMODE = 16)
;
; 10-Jul-03	gh 	Initial version.
; 11-Jul-03	gh	Restrict plotted points to those with roll_quality >=144 and > 2 blips
; 14-Jul-03 gh	Refine calculation of offset to account for variable period
; 23-Jul-03 gh	Change minimum quality for plot=1,4 mode from 144 to 216.
; 27-Jul-03 gh	Reorganize plotmode bits
;
npt = N_ELEMENTS(dbsolution)
t0 = dbsolution[0].sctime
max_offset = FLTARR(npt)
time = DBLARR(npt)
trel = dbsolution.sctime - t0
;
; PLOMODE=1 or 2 ==> AV_OFFSET vs time and/or AV_OFFSET histogram
;
IF (plotmode AND 3) NE 0  THEN BEGIN
	ok = WHERE(dbsolution.roll_quality GE 216 AND dbsolution.blipcount GT 2, nok)		; only plot the good points
	FOR n=0L,npt-1 DO BEGIN
		js = dbsolution[n].starcount
		max_offset[n] = MAX(ABS(dbsolution[n].stars[*].av_offset))*1000.		; converted to milliradians
		time[n] = HSI_SCTIME2ANY(dbsolution[n].sctime)
	ENDFOR
	max_offset 	= max_offset[ok]
;	time		= time[ok]
	log_max_offset = ALOG10( (max_offset < 1000.) > 0.01)
	IF (plotmode AND 1) NE 0 THEN $
		PLOT, trel[ok], (max_offset < 1000.) > 0.01, XTITLE='Relative Time (s)', YTITLE='AV_OFFSET (mrad)', /YLOG
	IF (plotmode AND 2) NE 0 THEN $
		PLOT, HISTOGRAM(log_max_offset, BINSIZE=0.1)>0.1,/YLOG, PSYM=6, XRANGE=[0,40]
ENDIF
;
; PLOTMODE=4 ==> Orbit-to-orbit phase shift vs time
;
IF (plotmode AND 4) NE 0 THEN BEGIN
	orbital_period = 96*60L
	norb = FIX(trel[npt-1]/orbital_period)						; number of complete orbits since t0
	expected_droll = DBLARR(norb)
	ok = WHERE(dbsolution.roll_quality GE 144)
	t = (INDGEN(norb)+0.3)*orbital_period						; Times (relative to t0) at which to interpolate solution
	roll = INTERPOL(dbsolution[ok].roll_phase, trel[ok], t)		; Returns roll angle (radians) once per orbit
	FOR n=0, norb-1 DO BEGIN
		i = WHERE(trel LE t[n]  AND trel GE t[n]-orbital_period $
								AND dbsolution.roll_quality GE 144, ni)		; Original indices between t
		IF ni EQ 0 THEN CONTINUE
		eff_period = MEAN (dbsolution[i].roll_period, /DOUBLE)
		expected_droll[n] = orbital_period / eff_period
	ENDFOR
	droll = ((roll - SHIFT(roll,1))/(2*!PI*1.E6)  - expected_droll+9999.5) MOD 1. - 0.5		; units of rotations
	PLOT, t/orbital_period, droll, psym=-6,   xrange=[0, 300]
	log_droll = ALOG10( ABS(droll)*2000.*!PI < 10000.)						; log(|droll|) in mrad
	PLOT, HISTOGRAM(log_droll, BINSIZE=0.1)>0.1, /YLOG, PSYM=6
ENDIF
;
; PLOTMODE=8 ==> Roll period vs time
;	PLotted squares are based on individual roll_period values.
;	Plotted x's are based on phase differences between successive 'ok' data points.
;
IF (plotmode AND 8) NE 0 THEN BEGIN
	ok = WHERE(dbsolution.roll_quality GE 216, nok)
	IF nok GE 2 THEN BEGIN
		time = trel[ok]
		rollperiod = dbsolution[ok].roll_period							; in seconds
		rotphase = dbsolution[ok].roll_phase * 0.5D-6 / !DPI			; in rotations
		xlabel = STRING(t0, FORMAT='("SC seconds since ", I11)')
		sctime0 = {seconds: t0, bmicro: 0L}
		timelabel = 'PMTRAS_DBASE_PLOT starting ' + HSI_SCTIME2ANY(sctime0, /STIME)
		PLOT, time, rollperiod, YSTYLE=16, PSYM=6, SYMSIZE=0.5, $
			XTITLE=xlabel, YTITLE='Roll Period(s)', TITLE=timelabel
		eff_time = (time + SHIFT(time, -1))/2
		dtime = SHIFT(time, -1) - time
		rotations = dtime *2.D / (rollperiod + SHIFT(rollperiod, -1)) - SHIFT(rotphase, -1) + rotphase
		nrot = ROUND(rotations)
		estperiod = dtime / (nrot - rotphase + SHIFT(rotphase, -1))
		OPLOT, eff_time[0:nok-2], estperiod[0:nok-2], PSYM=7, SYMSIZE=0.5
	ENDIF
ENDIF
;
; PLOTMODE=16 ==> Est_error vs time
;
IF (plotmode AND 16) NE 0 THEN BEGIN
	ok = WHERE(dbsolution.roll_quality GE 63)
	PLOT, trel[ok], dbsolution[ok].est_error*1000.,YSTYLE=16, XTITLE='Relative time(s)', YTITLE='Est_error(mr)'
ENDIF
RETURN
END