;+
;function: thm_interpolate_state
;
;Purpose: interpolates the low res STATE file
;
;         all variables are structures as produced by get_data
;
;keywords:
;
;
;Examples:
;      tha_spinper_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_spinper=thx_spinper) --> linear interpolation
;      tha_spinphase_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_spinper=thx_spinper,thx_spinphase=thx_spinphase) --> phase constructed according to the nearest neighbor spin phase, spin period
;      tha_spinras_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_spinras=thx_spinras) --> linear interpolation
;      tha_spindec_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_spindec=thx_spindec) --> linear interpolation
;      tha_spinalpha_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_spinalpha=thx_spinalpha) --> linear interpolation
;      tha_spinbeta_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_spinbeta=thx_spinbeta) --> linear interpolation
;      tha_pos_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_pos=thx_pos) --> spline interpolation
;      tha_vel_highres=thm_interpolate_state(thx_xxx_in=thx_xxx_in,thx_vel=thx_vel) --> spline interpolation
;
;Notes: under construction!!
;
;Written by Hannes Schwarzl
; $LastChangedBy: kenb-mac $
; $LastChangedDate: 2007-02-02 14:32:30 -0800 (Fri, 02 Feb 2007) $
; $LastChangedRevision: 278 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/thmsoc/trunk/idl/themis/state/thm_interpolate_state.pro $
;-


function thm_interpolate_state,thx_xxx_in=thx_xxx_in,thx_spinper=thx_spinper,thx_spinphase=thx_spinphase,thx_spinras=thx_spinras,thx_spindec=thx_spindec,thx_spinalpha=thx_spinalpha,thx_spinbeta=thx_spinbeta,$
		thx_pos=thx_pos,thx_vel=thx_vel


timeV=thx_xxx_in.X

is_kew_comb=0

if keyword_set(thx_xxx_in) && keyword_set(thx_spinper) then begin
	is_kew_comb=1
	;linearly interpolate the spinperiod
	sperInterp = interpol( thx_spinper.Y,thx_spinper.X,thx_xxx_in.X)

	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',sperInterp,'V',000.0)


endif

if keyword_set(thx_xxx_in) && keyword_set(thx_spinphase) && keyword_set(thx_spinper) then begin
	is_kew_comb=1
	;set a constant depending if the phase decreases with time
    ;phaseDir=-1.0d   ;decreasing phase
     phaseDir= 1.0d  ;increasing phase

	;interpolate the spin phase
	count=SIZE(thx_xxx_in.X,/N_ELEMENTS)


	countPhase=SIZE(thx_spinphase.X,/N_ELEMENTS)


	;create instant phase from closest phase information (nearest neighbor)
	; have not found a way to do nearest neighbor inperpolation using a built in IDL function??
	i=0
	timei=thx_spinphase.X[i]
	phase=dblarr(count)
	FOR j=0L,count-1L DO BEGIN
		; update header information if necessary (WHILE loop should work for all possible cases)
		WHILE ((thx_xxx_in.X[j] gt timei) && (i lt countPhase-1)) do begin
			i=i+1
			timei=thx_spinphase.X[i]
		ENDWHILE

		diffNext=thx_xxx_in.X[j]-thx_spinphase.X[i]
		IF (i gt 0L) THEN BEGIN
			diffPrev=thx_xxx_in.X[j]-thx_spinphase.X[i-1L]
			dnext=abs(diffNext)
			dprev=abs(diffPrev)
			IF (dnext lt dprev) THEN BEGIN
				closestPhase=i     ;next phase is closer
				diffAct=diffNext
			ENDIF ELSE BEGIN
				closestPhase=i-1L  ;previous phase is closer
				diffAct=diffPrev
			ENDELSE
		ENDIF ELSE BEGIN
			closestPhase=i ; no previous phase
			diffAct=diffNext
		ENDELSE


		;extract phase at time of interest
		phasei=thx_spinphase.Y[closestPhase]
		spineriodi=thx_spinper.Y[closestPhase]
		phaseRate=360.d0/spineriodi ;deg per sec
		phase[j]=(phasei+phaseRate*phaseDir*diffAct) mod 360.d0
		IF phase[j] lt 0.d0 THEN BEGIN
			phase[j]=phase[j]+360.d0
		ENDIF
	ENDFOR

	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',phase,'V',000.0)

endif



if keyword_set(thx_xxx_in) && keyword_set(thx_spinras) then begin
	is_kew_comb=1
	;linearly interpolate the right ascencion angle
	rasInterp = interpol( thx_spinras.Y,thx_spinras.X,thx_xxx_in.X)
	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',rasInterp,'V',000.0)
endif

if keyword_set(thx_xxx_in) && keyword_set(thx_spindec) then begin
	is_kew_comb=1
	;linearly interpolate the elevation angle
	decInterp = interpol( thx_spindec.Y,thx_spindec.X,thx_xxx_in.X)
	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',decInterp,'V',000.0)
endif



if keyword_set(thx_xxx_in) && keyword_set(thx_spinalpha) then begin
	is_kew_comb=1
	;linearly interpolate the alpha angle
	alphaInterp = interpol( thx_spinalpha.Y,thx_spinalpha.X,thx_xxx_in.X)
	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',alphaInterp,'V',000.0)
endif

if keyword_set(thx_xxx_in) && keyword_set(thx_spinbeta) then begin
	is_kew_comb=1
	;linearly interpolate the beta angle
	betaInterp = interpol( thx_spinbeta.Y,thx_spinbeta.X,thx_xxx_in.X)
	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',betaInterp,'V',000.0)
endif



if keyword_set(thx_xxx_in) && keyword_set(thx_vel) then begin
	is_kew_comb=1
	;spline interpolate the velocity
	count=n_elements(thx_xxx_in.X)
	velInterp=indgen(count,3,/float)
	velInterp(*,0) = interpol( thx_vel.Y,thx_vel.X,thx_xxx_in.X,/SPLINE)
	velInterp(*,1) = interpol( thx_vel.Y,thx_vel.X,thx_xxx_in.X,/SPLINE)
	velInterp(*,2) = interpol( thx_vel.Y,thx_vel.X,thx_xxx_in.X,/SPLINE)
	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',velInterp,'V',000.0)
endif

if keyword_set(thx_xxx_in) && keyword_set(thx_pos) then begin
	is_kew_comb=1
	;spline interpolate the position
	count=n_elements(thx_xxx_in.X)
	posInterp=indgen(count,3,/float)
	posInterp(*,0) = interpol( thx_pos.Y(*,0),thx_pos.X,thx_xxx_in.X,/SPLINE)
	posInterp(*,1) = interpol( thx_pos.Y(*,1),thx_pos.X,thx_xxx_in.X,/SPLINE)
	posInterp(*,2) = interpol( thx_pos.Y(*,2),thx_pos.X,thx_xxx_in.X,/SPLINE)
	thx_xxx_out = CREATE_STRUCT('X',timeV ,'Y',posInterp,'V',000.0)
endif

if	is_kew_comb eq 0 then begin
	PRINT, 'wrong combination of input arguments'
	stop
endif

RETURN, thx_xxx_out
end