pro gapsmth_sbsp,nsmoothr,ftime,wdelw,fitww

;  routine to smooth the thermal drift along the slit taking into account
;  the time gaps.  Time stamp associated with the frame in progress during the
;  gap has time value before the gap started, but data is associated with the
;  end of the gap.  This routine adjusts the time stamp before smoothing.

;  2008-02-22 - BWL;  modifications to not crash when it encounters a very
;		short segment of length <3


nfile = sizeof_sbsp(ftime,1)
fitww = fltarr(nfile)

igap = lonarr(2500)

; get median time in hours between the observations
;  take median time because there can be large gaps
;dtt = (ftime(nfile-1)-ftime(0))/float(nfile)
dtt = median(abs(ftime -shift(ftime,1)))

;  smoothing should be done on a fixed time interval (4 min or less),
;  not a fixed number of points.  The given value is typical of a 
;  normal map
nsmoot = fix(nsmoothr*4.8/(3600.*dtt))

;  go through again looking for discontinuities.  Accept as discontinuity
;  if time step jump > 8 average time between observations.  Record the 
;  index of the point IMMEDIATELY BEFORE the gap in time occurs
kntr = 0
for kk = 1,nfile-1 do begin
	if (abs(ftime(kk) - ftime(kk-1))) gt 8.*dtt  then begin
		igap(kntr) = kk-1
		kntr = kntr + 1
	endif
endfor
;  make gap counter only as long as number of gaps
if kntr gt 0 then begin
	igap = igap(0:kntr-1)
endif else begin
	igap = 0
endelse
ngap = kntr
	
;  smooth only up through time 2 steps before the apparent beginning
;  of time gap since this gap is usually caused by a hold of the SP program,
;  and one measurement is started, but not finished until the end of the gap
istrt=0
if ngap gt 0 then begin
	for kgap = 0,ngap-1 do begin
		iend = igap(kgap) - 2
;  check for short segments
		if iend lt istrt then iend = istrt
		nsegmt = iend-istrt+1
		;nsmthw = min([nsmoothr,nsegmt/2])
;  smoothing on a fixed interval of time
		nsmthw = min([nsmoot,nsegmt/2])
		temp = wdelw(istrt:iend)
		fitww(istrt:iend) = smooth(temp,nsmthw,/edge_truncate)
;  just duplicate the smoothed curve at the end of the segment
;  if short segment, then duplicate only to the end of the segment
		iend2 = iend+2
		ibegn = iend+1
		if iend2 gt igap(kgap) then iend2 = igap(kgap)
		if ibegn gt iend2 then ibegn = iend2
		fitww(ibegn:iend2) = fitww(iend)
;  adjust time for observation point immediately adjacent to time gap
		if kgap gt 0 then ftime(istrt) = ftime(istrt+1) - dtt
;  increment the start index
		istrt = igap(kgap)+1
	endfor

;  smooth the last segment
	iend = nfile-1
; if only one gap, adjust the first point observation time
	if ngap eq 1 then ftime(istrt) = ftime(istrt+1) - dtt
	nsegmt = iend-istrt+1
	temp = wdelw(istrt:iend)
	nsmthw = min([nsmoot,nsegmt/2])
	fitww(istrt:iend) = smooth(temp,nsmthw,/edge_truncate)
endif else begin
;  case of no gaps
;  avoid endpoints by duplicating adjacent points before smoothing
	nsegmt = nfile
	if nsegmt ge 2 then begin
		nsmthw = min([nsmoot,nsegmt/2])
		tempww = wdelw
		tempww(0) = tempww(1)
		tempww(nfile-1) = tempww(nfile-2)
	endif else begin
		tempww = wdelw
	endelse
	fitww = smooth(tempww,nsmthw,/edge_truncate)
endelse



if nsmthw le 1 then goto,jump22
;  condition ends of segments -- smooth does not do a good job
;  if end points of segments differ substantially from mean

;  fix is to cycle through segments, conditioning ends with 
;  smooth transition to a linear fit over the lesser of the 
;  smoothing parameter or the number of points in the segment

;  first find the start, end points of the segments
;  number of segments = ngap+1
strtseg = lonarr(ngap+1)
endseg = lonarr(ngap+1)
strtseg(0) = 0
if ngap gt 0 then begin
	for kgap = 0,ngap-1 do begin
		endseg(kgap) = igap(kgap)
		strtseg(kgap+1) = igap(kgap) + 1
	endfor
endif
endseg(ngap) = nfile-1

;  condition ends of segments
for kgap = 0,ngap do begin
	istrts = strtseg(kgap)
        iends = endseg(kgap)
	nsegmt = iends - istrts + 1
	nlinr = min([nsmoot,nsegmt])
	nl1 = nlinr-1
;  for case of less than 3 elements, bypass
	if nsegmt le 2 then goto,jump11
;  must have two points for linear fit
	if nl1 le 1 then goto,jump11

;  linear trend at start of segment, dont fit first element
	xxx = ftime(istrts+1:istrts+nl1) 
	yyy = wdelw(istrts+1:istrts+nl1)
	coef = linfit(xxx,yyy,yfit=lfit)
;  skip fitting first point, extrapolate to that
	lff = fltarr(nlinr)
	lff(1:*) = lfit
	lff(0) = lfit(0) - (ftime(istrts+1)-ftime(istrts))*coef(1)
	lfit = lff
	for ii = 0,nl1 do begin
		blendfac = (ftime(istrts+nl1)-ftime(istrts+ii))/ $
				(ftime(istrts+nl1)-ftime(istrts))
		fitww(istrts+ii) = blendfac*lfit(ii) + $
			(1. - blendfac)*fitww(istrts+ii)
	endfor

;  linear trend at end of segment
	xxx = ftime(iends-nl1:iends-1)
	yyy = wdelw(iends-nl1:iends-1)
	coef = linfit(xxx,yyy,yfit=lfit)
;  skip fitting last point, extrapolate to that
	lff = fltarr(nlinr)
	lff(0:nl1-1) = lfit
	lff(nl1) = lfit(nl1-1) + (ftime(iends)-ftime(iends-1))*coef(1)
	lfit = lff
	for ii = 0,nl1 do begin
		blendfac = (ftime(iends-nl1+ii)-ftime(iends-nl1))/ $
			(ftime(iends)-ftime(iends-nl1))
		fitww(iends-nl1+ii) = blendfac*lfit(ii) + $
			 (1.-blendfac)*fitww(iends-nl1+ii)
	endfor
	jump11:
endfor
jump22:


return
end