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