PRO	PMTRAS_DBASE_INTERPOLATE, dbsolution, ndeg, gapcounts, return_flag
;
; Interpolates missing or marginal points in dbsolution with a polynomial of degree, ndeg.
; Sets return_flag=1 if any changes to dbsolution were made; 0 otherwise.
; Loads gapcounts array with gap statistics.
;
; 21-Aug-03 gh	Initial version.
; 30-Nov-03 gh	Revise formula for calculation of est_error.
;				Ease fit tolerance from 1 to 1.3 arcminutes (per gh notebook analysis 11-30-03).
;
return_flag 	= 0
npt 			= N_ELEMENTS(dbsolution)
fit_tolerance 	= 0.0004								; Fit tolerance = 1.3 arcminutes
twopi 			= 2.*!DPI
;
; Identify good data points in original dbsolution.
igood		= WHERE(dbsolution.roll_quality GE 192, ngood)
IF (ngood LT 2) THEN BEGIN
	PRINT, 'PMTRAS_DBASE_INTERPOLATE: Fails with only ', ngood, 'data points.'
	RETURN
ENDIF
scgoodtime	= dbsolution[igood].sctime
gapcounts	= INTARR(4)					; gapcounts[0] = total number of gaps
;										; gapcounts[1] = number of short gaps (< 800s)
;										; gapcounts[2] = number of short gaps with good coverage
;										; gapcounts[3] = number of successfully filled gaps.
;
; Identify prospective gaps.
dt				= SHIFT(scgoodtime, -1) - scgoodtime	; good array index corresponding to earlier time, dt >0
dummy			= WHERE (dt GT 65, ng)					; only used for getting a gapcount
jgap			= WHERE (dt GT 65 AND dt LT 800, ngap)	; selects indices in good array of start of 65-800s gaps
gapcounts[0]	= ng
gapcounts[1]	= ngap
IF (ngap EQ 0) THEN BEGIN
	PRINT, 'PMTRAS_DBASE_INTERPOLATE: No short gaps to fill.'
	RETURN
ENDIF
igapstart		= igood[jgap]								; original array indices at start of gaps
igapend			= igood[jgap+1]								; original array indices at end of gaps
missing			= dt[jgap]/64-1								; array of number of missing data points
;
; Define a structure for a single interpolated data point.
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
dbinterpoint = {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:	192B,	$	; 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	}
;
; Begin loop over individual gaps.
; First task is to select the good 'reference' points that will be used to fit solution vs time.
; Points on each side of gap should be the minumum satisfying the 3 criteria that:
;	A. There must be at least ndeg+1 points;
;	B. There must be at least a number equal to the number of missing points in the gap
;	C. The separation of at least 1 of the points from the nearest edge of the gap be at least as large as half the
;			width of the gap.
FOR n=0, ngap-1 DO BEGIN
	jgapnstart 	= jgap[n]									; index in good array of start of this gap
	jgapnend	= jgap[n]+1									; index in good array of end   of this gap
	missingn	= missing[n]								; number of missing data points in this gap
	dtgapn		= dt[jgap[n]]
	nmin		= missingn > (ndeg + 1)						; minimum number of 'before' points reqd
; Identify relevant pregap indices
	jp0			= jgapnstart - INDGEN(nmin)					; indices in good array of points shortly before this gap.
	IF (jp0[nmin-1] LT 0 ) THEN CONTINUE					; insufficient points to interpolate - go to next gap
	dtp0	= scgoodtime[jgapnstart] - scgoodtime[jp0]		; array of (increasing +ve) intervals until start of gap
	m0		= WHERE(dtp0 GE dtgapn/2.)						; points far from gap
	mpre	= (ndeg+1) > (m0[0]+1) > missingn				; final number of 'before' points
	jpre	= INDGEN(mpre) + jgapnstart - mpre + 1			; indices of relevant pregap points in good array
; Identify relevant postgap indices
	jp1		= INDGEN(nmin) + jgapnend						; indices in good array of points shortly after this gap.
	IF (jp1[nmin-1] GT ngood) THEN CONTINUE 				; insufficient points to interpolate - go to next gap
	dtp1	= scgoodtime[jp1] - scgoodtime[jgapnend]		; array of (increasing +ve) intervals until start of gap
	m1		= WHERE(dtp1 GE dtgapn/2.)						; points far from gap
	mpost	= (ndeg+1) > (m1[0]+1) > missingn				; final number of 'before' points
	jpost	= INDGEN(mpost) + jgapnend						; indices of relevant postgap points in good array
; Combine to specify reference points for this gap.
	jref	= [jpre, jpost]									; array of indices in good array of reference points
	mref	= mpre + mpost									; number of reference points (>=8).
	x = scgoodtime[jref]									; = times in seconds
	IF x[mref-1] - x[0] GT 2400 THEN CONTINUE				; range of reference times is too great - go to next gap
	gapcounts[2] = gapcounts[2] + 1
;
; Polynomial fit of points around the gap.
	soln = hsi_pmtras_dbase_expand (dbsolution[igood[jref]])	; returns normal roll_solution structure
	t0f	= soln.t0_ref.seconds									; sctime.seconds of first datum
	x = soln.reltime
	y = soln.posn_angle											; radians
	fitparm = POLY_FIT(x, y, ndeg, /DOUBLE, YFIT=yfit)			; polynomial fit
; 			est_error = SQRT(TOTAL((y-yfit)^2)/(mref-ndeg-1))	; Obsolete formula replaced 11-30-03 gh.
 	est_error = SQRT( TOTAL((y-yfit)^2) / (mref*(mref-ndeg-1)))
	IF (est_error GT fit_tolerance) THEN CONTINUE				; scatter is too large - go to next gap
	gapcounts[3] = gapcounts[3] + 1
;
; Flag existing gap points in dbsolution for deletion.
	igapstartn = igapstart[n]
	igapendn	= igapend[n]
	FOR i = igapstartn+1, igapendn-1 DO dbsolution[i].sctime = -1L
;
; Calculate solutions for all missing points.
	t0m = dbsolution[igapstartn].sctime +64 - t0f				; time of 1st missing point in same units as 'x'
	t = t0m + 64L*INDGEN(missingn)								; times of all the missing points.
	ygap = DBLARR(missingn)
    FOR k=0, ndeg DO ygap = ygap + fitparm[k]*t^k	  		    ; array of interpolated values (increasing radians)
	roll_phase 	= LONG((ygap MOD twopi)*1.D6)	; microradians
	ygap0								= [y[mpre-1], ygap]		; add last pregap point to interpolated value array
	xgap0								= [x[mpre-1], t]
	roll_period		= twopi * (xgap0[1:missingn] - xgap0[0:missingn-1]) / $
												  (ygap0[1:missingn] - ygap0[0:missingn-1])
	approx_angvel		= BYTE(128 * twopi / roll_period)
;
; Load a structure with one point at a time and append it to an array of new points.
	FOR m=0, missingn-1 DO BEGIN
		dbinterpoint.sctime			= t[m] +t0f
		dbinterpoint.roll_phase 	= roll_phase[m]
		dbinterpoint.approx_angvel 	= approx_angvel[m]
		dbinterpoint.roll_period 	= roll_period[m]
		dbinterpoint.est_error 		= est_error
		IF N_ELEMENTS(dbnew) EQ 0 THEN dbnew = dbinterpoint ELSE dbnew = [dbnew, dbinterpoint]
	ENDFOR
ENDFOR
;
; Eliminate obsolute points in dbsolution, insert new points at their proper location.
IF N_ELEMENTS(dbnew) EQ 0 THEN RETURN						; no changes to database
iok	= WHERE(dbsolution.sctime GE 0)
dbtemp = [dbsolution[iok], dbnew]
iseq		= SORT(dbtemp.sctime)
dbsolution	= dbtemp[iseq]									; dbsolution is time-ordered again.
return_flag = -1
RETURN
END