PRO xrt_dustgrow,index,dust_out,dust,hist=hist, quiet=quiet
; =========================================================================
;+
;
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;
;       XRT_DUSTGROW 
;
; CATEGORY:
;
;       Hinode XRT data analysis.
;
; PURPOSE:
;
;       Makes an approximate estimate of the time dependent growth of 
;       semi-opaque material around dust specks.  
;
; CALLING SEQUENCE:
; 
; 
;      XRT_DUSTGROW,index,dust_out,[dust],[hist=hist]
;
; INPUTS:
;
;      INDEX        - [Mandatory] (a standard XRT index format)
;                      Index structure for the data for which dust
;                      growth estimate is requested.
;
;      DUST         - [Optional] full frame 1x1 dust image 
;                      (otherwise program will get it from dust_pix1.dat)
;
;
;   OUTPUTS:
;
;       DUST_OUT    - [Mandatory] Full frame frame 1x1 dust image corresponding
;                      to input index, modified for dust growth (if any) 
;
;  KEYWORDS: 
;       HIST        - [Optional] output history addednum note (can be added to
;                      index.history if user desires) 
;      QUIET        - run without printing [default - not quiet]

;
;
; EXAMPLES:
;
;       Basic call to get modified dust image 
;        IDL> xrt_dustgrow, index, dust_out
;
;       Call to get a modified dust image, using input dust image and
;         getting history update note 
;      IDL> xrt_dustgrow, index, dust_out, dust, hist=hist_grow
;
;  
; COMMON BLOCKS:
;
;       None
;
; NOTES:
;
; Uses a calibration Dave McKenzie's measurements of growth of the largest
; dust grain to dilate the original dust map suitably.  Sets a (large) constant
; correction if outside calibrated time range and prints a warning.  
; The model for the "dust growth" will likely be improved over time.
;
;
; MODIFICATION HISTORY:
;
;      31-Oct-2013 - S.Saar (written)
;      20-Nov-2013 - S.Saar (fixed path to darks if dark not predefined)
;      03-Feb-2014 - S.Saar (changed path to gband data file, added warning
;                           if past calibrated period, fixed constant [large]
;                            size for growth in this regime, changed note above)
;      14-Mar-2014 - S.Saar (commented out path used for local testing)
;
;-
; =========================================================================

progver = 'v2013-Oct-31' ;--- (S. Saar)
progver = 'v2013-Nov-20' ;--- (S. Saar)
progver = 'v2014-Feb-03' ;--- (S. Saar)
progver = 'v2014-Mar-14' ;--- (S. Saar)

q_qt=keyword_set(quiet)

date=index.date_obs
tmx=max(anytim(date))
path0=getenv('SSW_XRT')+'/idl/util/'
path1=getenv('SSW_XRT')+'/idl/response/contam/'

ndust=n_elements(dust)
if ndust eq 0 then begin
   dfile='dust_pix1.dat'             ; read in dust file
   dlen=1035
   dust =intarr(2048,2048)
   openr,1,path0+dfile
   arr=intarr(2,dlen)
   readf,1,arr
   close,1
   xx=reform(arr(0,*))                ; read in x & y coords of features
   yy=reform(arr(1,*))
   dust(xx,yy)=1                       ; construct integer full-frame image
endif

glarge = 200.                               ; gband estimate after last 
                                            ; calibration point 
tdustgrow='2012-08-10T00:00:00.00Z'         ; dust growth starts after eclipse
                                            ; season after light leak
tdgrow=anytim(tdustgrow)
; path1=''                                  ; (for local testing)
kern=[[0,1,0],[1,1,1],[0,1,0]]              ; growth kernel
dustg = dilate(dust,kern)                  
dustg = dilate(dustg,kern)                  ; dilate dust map twice as a minmum 
if tmx ge tdgrow then begin                 ; if in dust growth epoch...
   restore,path1+'gband_vs_bakeout.sav'     ; restore gband time variation
   tgband=[tbegin,tend]
   tgband=anytim(tgband)
   tdbmax=max(tgband)                       ; max calibrated date
   if tmx le tdbmax then begin           ; in calibrated range? find gband level
      ist=sort(tgband)
      tgband=tgband[ist]
      gband=[gband_begin,gband_end]
      gband=gband[ist]
      gbnd=interpol(gband,tgband,tmx)       ; interpolate Gband level
   endif else gbnd=glarge                   ; or set to large constant
   cws=[-1.019,0.0238]
   wspot=round(poly(gbnd,cws))          ; use width_spot(Gband level) fit
   for j=1,wspot do dustg = dilate(dustg,kern)             ; dilate dust map
   if tmx le tdbmax then begin                           ; in calibrated range? 
      if wspot gt 0 then hist='Dust growth estimated ' + $    ; if grown & in 
        'via dust dilation by '+strtrim(string(fix(wspot)),2)+ $ ; calibrated 
        ' pixels  (1x1).' else hist=''              ; time range, add to history
   endif else begin                                      ; if outside cal times 
      hist='Warning: outside calibrated time range for dust growth. Dust map'+ $
       ' arbitrarily corrected by '+strtrim(string(fix(wspot)),2)+' pixels ' + $
       '(1x1).'                                         ; add warning to history
      if q_qt ne 1 then begin                            ; and print warning 
         print,'Warning: outside calibrated time range for dust growth. '+ $
          'Dust map arbitrarily dilated by '+strtrim(string(fix(wspot)),2)+ $
          ' pixels  (1x1).' 
      endif
   endelse
endif else hist=''                     ; or if no growth, leave blank

dust_out=dustg
return
end
;
;
;