;+
; $Id: resamplemap.pro,v 1.1 2010/03/23 16:17:44 nathan Exp $
;
; PROJECT:
;  SOHO-LASCO
;
; NAME:
;  RESAMPLESMAP
;
; PURPOSE:
;  Resample a synoptic or slice map on a regularly gridded time 
;
; CATEGORY:
;  global activity, interpolation
;
; DESCRIPTION: 
;
; INPUTS:
;  map : 2D array, the time is the X axis
;  time : time corresponding to the input array
;  regtime : the time array where the user wants to interpolate the map
;
; OUTPUTS:
;  mapbinned : re-sampled map obtain with the binned method. The size
;              of the bin is equal to to size of a pixel.
;  mapreg : re-sampled map obtain with the IDL spline function. Too
;           long gaps are filled with value 0 according to
;           halfinterpsize keyword.
;  mapreginterp : re-sampled map obtain with the IDL spline function.
;
; INPUT KEYWORDS:
;  halfinterpsize : half width for the interpolation, in time
;                   unit. Outside this boundary the signal is filled with 0. 
;
; OUTPUT KEYWORDS:
;
; CVSLOG:
; $Log: resamplemap.pro,v $
; Revision 1.1  2010/03/23 16:17:44  nathan
; moved from dev/synomap
;
; from resamplemap1.pro,v
; 
; Revision 1.1  2007/01/24 21:57:50  thernis
; First commit
;
;-
pro resamplemap,map,time,regtime,mapbinned,mapreg,mapreginterp,halfinterpsize=halfinterpsize

;time=lst.jdint-lst(0).jdint+lst.jdfrac-lst(0).jdfrac
;dobs2jd,'1996/09/01','00:00:00',jdint1,jdfrac1
;dobs2jd,'1996/09/30','23:59:59',jdint2,jdfrac2
;regtime=findgen(1000)/999*(jdint2-jdint1+jdfrac2-jdfrac1)+((jdint1-lst(0).jdint)+(jdfrac1-lst(0).jdfrac))


if n_elements(halfinterpsize) eq 0 then halfinterpsize=0.5 ; -- max gap in days for interpolation


sx=(size(map))(1)
sy=(size(map))(2)


stimeout=n_elements(regtime)

mapreg=fltarr(stimeout,sy)
mapbinned=fltarr(stimeout,sy)


; --- re-sample by interpolation of each row of the map
for i=0,sy-1 do begin
    r=interpol(map(*,i),time,regtime)
    mapreg(*,i)=r
endfor

; --- re-sample by the bin method
; --- and clean up the fake data from the interpolation
timederiv=time-shift(time,1)
mg=where(timederiv gt halfinterpsize*2,cntg)
mapreginterp=mapreg
halfbinsize=.5*(regtime(stimeout-1)-regtime(0))/float(stimeout-1) ; -- day per pix
for i=0,stimeout-1 do begin
    diff=time-regtime(i)
    mm=where(diff le halfbinsize and diff ge -halfbinsize,cnt)
    if cnt gt 0 then mapbinned(i,*)=total(map(mm,*),1)/cnt 
    mm=where(diff le halfinterpsize and diff ge -halfinterpsize,cnt)
    if cnt eq 0 then mapreg(i,*)=0

endfor

; --- cut at the hedges of the hole of data to be cleaner
sxmap=(size(mapreg))(1)
for i=0,cntg-1 do begin
    m1=where(regtime le time(mg(i)-1),cnt1)
    if cnt1 gt 0 then m1=m1(cnt1-1)
    m2=where(regtime ge time(mg(i)),cnt2)
    if cnt2 gt 0 then m2=m2(0)
    if cnt1 gt 0 and cnt2 gt 0 then mapreg(m1 > 0 :m2 < (sxmap-1),*)=0

endfor

return
end