PRO PMTRAS_DBASE_FIT_FUNCTION, x, a, f, pder
;
; Special purpose procedure used with CURVEFIT to evaluate
; f(x) = A[0] + A[1]*x + A[2]*x^2 + A[3]*COS(2PIx/86400) + A[4]*SIN(2PIx/86400) and its derivatives
;
twopiday = 2*!DPI / 86400.
fourpiday = 2*twopiday
npt = N_ELEMENTS(x)
f 		= a[0] + a[1]*x + a[2]*x^2 + a[3]*COS(twopiday*x)  + a[4]*SIN(twopiday*x) 	$
								   + a[5]*COS(fourpiday*x) + a[6]*SIN(fourpiday*x)
pder 	= [[FLTARR(npt) + 1.], [x], [x^2], [COS(twopiday*x)],   [SIN(twopiday*x)],	$
										   [COS(fourpiday*x)],  [SIN(fourpiday*x)]]
END
PRO PMTRAS_DBASE_FIT, dbsolution, fitparm, INFLAG=inflag
;
; Does a multiparameter fit to roll data in structure array with the master database format.
;
; /INFLAG indicates use of the alternate input format as generated by hsi_pmtras_lookup.
;
; Time-dependence of roll angle is assumed to be of the form:
;	y(t) = K + At + Bt^2									; offset, roll rate, angular acceleration
;			+ C1*cos(2pi*t  /P) + S1*sin(2pi*t/P)			; orbital period
;           + C2*cos(2pi*t*2/P) + S2*sin(2pi*t*2/P)			; 2nd harmonic of orbital period
;           + C3*cos(2pi*t*2/P) + S3*sin(2pi*t*2/P)			; 3rd harmonic of orbital period
;           + C4*cos(2pi*t*2/P) + S4*sin(2pi*t*2/P)			; 4th harmonic of orbital period
;           + C0*cos(2pi*t* /D) + S0*sin(2pi*t* /D)			; diurnal variation
;
; Method: After an initial removal of polynomial term, the orbital period will be determined by fitting the orbital term.
;		  This orbital period will then be treated as a constant.
;		  Next, the harmonics of the orbital periods will be fit, with their parameters then treated as constants.
;		  Next, the orbital terms will be subtracted from the original data.
;		  Finally the polynomial terms and the diurnal variation will be fit using CURVEFIT.
;
; Output units are degrees and seconds (relative to t0_ref).
;
; 14-Apr-04	gh 	Initial version.
; 13-Oct-05 gh	Add /INFLAG keyword

;
IF KEYWORD_SET(inflag) EQ 0 THEN BEGIN
	ok 			= WHERE(dbsolution.roll_quality GE 192)
	solnok 		= hsi_pmtras_dbase_expand(dbsolution[ok])	; convert to structure with monatonically increasing radians
	t0 			= dbsolution[0].sctime				; LONG
	trel 		= dbsolution.sctime - t0			; LONG array
	trelok 		= DOUBLE(trel[ok])
	phaseok		= solnok.posn_angle *!RADEG		; convert to degrees
ENDIF ELSE BEGIN
	t0 			= dbsolution.t0_ref.seconds			; LONG
	trelok		= dbsolution.reltime
	phaseok		= dbsolution.posn_angle *!RADEG
ENDELSE
;
; Initial polynomial fit.
fitparm1 	= POLY_FIT(trelok, phaseok, 2, YFIT = phasefit1, /DOUBLE)
resideg1	= phaseok - phasefit1 			; residuals after cubic 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)
IF a2 GE (a1 > a3) THEN orbperiod = orbperiod*(1 + deltap * 0.5*(a3-a1)/(2*a2-a3-a1))
;
; Fit terms that have period of 1 orbit and its harmonics
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)
;
; Do the final fit
twopiday	= 2*!DPI/86400.			; angular velocity with 1-day period
fitparm2	= [0., 0., 0., 0., 0., 0., 0.]				; starting values
yfit		= CURVEFIT(trelok, resideg5, FLTARR(N_ELEMENTS(trelok))+1., fitparm2, $
				FUNCTION_NAME='pmtras_dbase_fit_function', YERROR=yerror, /DOUBLE)
print, fitparm1
resid	= resideg5 - yfit
!P.MULTI= [0,1,2]				; 2 rows of plots
UTPLOT, trelok, resideg1, hsi_sctime2any({seconds: t0, bmicro:0}), YTITLE='Phase residual (degrees)', $
			PSYM=1, SYMSIZE=0.25, TITLE='ROLL PHASE RESIDUAL WITH QUADRATIC FIT'
UTPLOT, trelok, resid, hsi_sctime2any({seconds: t0, bmicro:0}), YTITLE='Phase residual (degrees)', $
			PSYM=1, SYMSIZE=0.25, TITLE='FINAL ROLL PHASE RESIDUAL'
avperiod 	= 360./(fitparm1[1]+fitparm2[1])
dperiod		= -avperiod^2/(360.) * (fitparm1[2]+fitparm2[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')"
fitparm 	= [fitparm1[0]+fitparm2[0], fitparm1[1]+fitparm2[1], fitparm1[2]+fitparm2[2], orbperiod, $ ; final fit parms
				orbparm1[2], orbparm1[3], orbparm2[2], orbparm2[3],								$
				orbparm3[2], orbparm3[3], orbparm4[2], orbparm4[3],								$
				fitparm2[3], fitparm2[4], fitparm2[5], fitparm2[6] ]
RETURN
END