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)
;	Bit 5 = 1 ==> Phase residual vs time				(eg PLOTMODE = 32)
;
; 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
; 16-Mar-04 gh	Implement Phase residual plot.
; 17-Mar-04 gh	Change minimum quality for plot 32 to 192 to accomodate RAS solutions.
; 12-Apr-04 gh	Adapt plot 32 code to new version of hsi_remove_sinusoid
;				Settle on order 3,6 polynomial fits for plot 32.
; 13-Apr-04 gh	Minor change to plot format.
;				Change to 2,6 polynomial fit for plot 32.
;				Add printed display of period and dperiod/dt for plot 32.
;
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 = dbsolution[ok].sctime - dbsolution[0].sctime
		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
	SAVE, estperiod[0:nok-2]
	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
;
; PLOTMODE=32 ==> Phase residual vs time
;
IF (plotmode AND 32 ) NE 0 THEN BEGIN
	ok 			= WHERE(dbsolution.roll_quality GE 192)
	solnok 		= hsi_pmtras_dbase_expand(dbsolution[ok])	; convert to structure with monatonically increasing radians
	trelok 		= trel[ok]
	phaseok		= solnok.posn_angle
	fitparm1 	= POLY_FIT(trelok, phaseok, 2, YFIT = phasefit1, /DOUBLE)
	resideg1	= (phaseok - phasefit1) * !RADEG			; residuals after cubic fit
	UTPLOT, trelok, resideg1, hsi_sctime2any({seconds: t0, bmicro:0}), PSYM=1, SYMSIZE=0.25, $
			YTITLE='Phase residual (degrees)', TITLE='Roll phase residual after quadratic fit'
;
;Find best fit period by fitting nominal period and +- 1%, and fitting amplitude.
	orbperiod	= 96*60
	deltap		= 0.01
	orbparm1	= hsi_remove_sinusoid(trelok, resideg1, orbperiod*(1-deltap), YRESID=test1)
	orbparm2	= hsi_remove_sinusoid(trelok, resideg1, orbperiod, 			  YRESID=test2)
	orbparm3	= hsi_remove_sinusoid(trelok, resideg1, orbperiod*(1+deltap), YRESID=test3)
	a1			= SQRT(orbparm1[2]^2 + orbparm1[3]^2)
	a2			= SQRT(orbparm2[2]^2 + orbparm2[3]^2)
	a3			= SQRT(orbparm3[2]^2 + orbparm3[3]^2)
	enable_period_fit = 1
	IF a2 GE (a1 > a3) AND enable_period_fit EQ 1 THEN orbperiod = orbperiod*(1 + deltap * 0.5*(a3-a1)/(2*a2-a3-a1))
	orbparm1		= hsi_remove_sinusoid(trelok, resideg1, orbperiod,   YRESID=resideg2, /DOUBLE)
	orbparm2		= hsi_remove_sinusoid(trelok, resideg2, orbperiod/2, YRESID=resideg3, /DOUBLE)
	orbparm3		= hsi_remove_sinusoid(trelok, resideg3, orbperiod/3, YRESID=resideg4, /DOUBLE)
	orbparm4		= hsi_remove_sinusoid(trelok, resideg4, orbperiod/4, YRESID=resideg5, /DOUBLE)
	fitparm2 		= POLY_FIT(trelok, resideg5, 6, YFIT = phasefit, /DOUBLE, YERROR=yerror)
	resideg			= resideg5 - phasefit
	UTPLOT, trelok, resideg, hsi_sctime2any({seconds: t0, bmicro:0}), YTITLE='Phase residual (degrees)', $
				PSYM=1, SYMSIZE=0.25, TITLE='ROLL PHASE RESIDUAL'
	avperiod 	= 2*!DPI/fitparm1[1]
	dperiod		= -avperiod^2/(2*!DPI) * fitparm1[2] * 96*60
;
; Print some parameters
	PRINT, avperiod, 		FORMAT="('Average spin period:         ', F9.6, ' seconds')"
	PRINT, dperiod,  		FORMAT="('Spin period change per orbit:', F9.6, ' seconds')"
	PRINT, orbperiod/60, 	FORMAT="('Orbital period estimate:     ', F9.5, ' minutes')"
	PRINT, yerror,			FORMAT="('RMS residual:                ', F9.3, ' degrees')"
ENDIF
RETURN
END