;+
; NAME:
;      test_stacking_ias_library
;
; PURPOSE:
;      tests and checks the installed IDL stacking IAS library
;
; CALLING SEQUENCE:
;      test_stacking_ias_library
;
; INPUTS:
;     none, except files on disk
;
; OUTPUTS:
;    none, except screen messages
;
; EXAMPLE:
;      test_stacking_ias_library
;      test_stacking_ias_library, /healpix
;
; MODIFICATION HISTORY:
;    07-Oct-2009: written by H. Dole and M. Bethermin, IAS
;    9-Dec-2009: Add test of fullstack_xy by M. Bethermin
;
;
;
;     Copyright Insitut d'Astrophysique Spatiale
;
;     This program is free software: you can redistribute it and/or modify
;     it under the terms of the GNU General Public License as published by
;     the Free Software Foundation, either version 2 of the License, or
;     (at your option) any later version.
;
;     This program is distributed in the hope that it will be useful,
;     but WITHOUT ANY WARRANTY; without even the implied warranty of
;     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;     GNU General Public License for more details.
;
;     You should have received a copy of the GNU General Public License
;     along with this program.  If not, see <http://www.gnu.org/licenses/>.
;
;-


PRO test_stacking_ias_library, healpix=healpix

;------------------------------------------
; XY part of the Stacking IAS Library
;------------------------------------------

;Generate a map to test fullstack_xy
  Nmap = 100
  Nmasked = 3000
  Ns = 5000
  map = fltarr(Nmap,Nmap) ; create map
  psf = [[0.,0.15,0.],[0.15,0.4,0.15],[0.,0.15,0.]] ; create very simple PSF normalized to unity
  x = floor(Nmap*randomu(seed,Ns)) ; random x position of INPUT sources
  y = floor(Nmap*randomu(seed,Ns)) ; random y position of INPUT sources
  bad = where(x eq Ns) ; make sure sources are within map
  if bad(0) ne -1 then x(bad) = Ns-1
  bad = where(y eq Ns)
  if bad(0) ne -1 then y(bad) = Ns-1
  for k = 0L, Ns-1 do begin
     map[x(k),y(k)] = map[x(k),y(k)]+1.
  endfor
  map = convolve(map,psf) ; convolution by PSF
  masked = floor((Nmap*1.*Nmap-1.)*randomu(seed,Nmasked))
  mask = map*0.+1.
  mask(masked) = 0.

;Test fullstack_xy
  halfsize = 15. ; image stack size
  sel = where(x gt halfsize and y gt halfsize and x lt Nmap-halfsize and y lt Nmap-halfsize)
  fullstack_xy, map, x(sel), y(sel), flux, size = halfsize, img_stack=stacktest, $
                psf = psf, unc_flux = unc_flux,Nboot=20, /silent, weight = mask
  
  print, '###########################################'
  if abs(flux-1.) lt 5.*unc_flux then print, 'fullstack_xy: OK' else print, 'fullstack_xy: FAILED'
  print, '###########################################'

factor = 5
LOADCT, 3
WINDOW, 1, xs=Nmap*factor, ys=Nmap*factor, title='Simulated Input Map'
TVSCL, REBIN(map, Nmap*factor, Nmap*factor, /sa)

WINDOW, 2, xs=(2*halfsize+1)*factor, ys=(2*halfsize+1)*factor, title='Stack Map'
TVSCL, REBIN(stacktest, (2*halfsize+1)*factor, (2*halfsize+1)*factor, /sa)


;------------------------------------------
; HEALPIX part of the Stacking IAS Library
;------------------------------------------
IF KEYWORD_SET(healpix) THEN BEGIN
  pathwmap = './check/'
  wmapfile = 'wmap_band_imap_r9_5yr_W_v3_temperature.save'
  RESTORE, pathwmap+wmapfile
;mollview, mapt, /online, /nest
  skyinput = mapt
  nest=1
; weight map equals to 1
  weight = skyinput*0. + 1.

; Output Stack pixel size and map size
; => CHANGE HERE
  reso_arcmin = 10.             ; in arcmin
  size = 20

; Lat, Long (in whatever system: C, G or E) coordinates of input catalog
; => CHANGE HERE catalog file
;  READCOL, pathwmap + 'clusters-z0-0.3T4-15_2.txt', ra, dec, F='F,F'
  READCOL, pathwmap + 'wmap_ptsrc_catalog_p3_5yr_v3.txt', ra, dec, l, b, K, Ka, Q, V, W, eK, eKa, eeQ, eV, eW ;, F='F,F, F, F, F, F, F, F, F, F, F, F, '

; select undetected sources in band W
  good = WHERE(2.*eW GE W, ng)
  PRINT, ' Number of selected sources in WMAP band W: ', ng
  l = l(good)
  b = b(good)

; Convert input catalog in the same system as input map
; here RA,DE->l,b

  phi = l/!RADEG
  theta = (90.-b)/!RADEG
  ANG2PIX_NEST,  512, theta, phi, pixels

  LOADCT, 5

; Stack
  DOSTACK_HEALPIX, skyinput, l, b, reso_arcmin, img_stack_output, size=size, datacube=datacube, Nstacked=Nstacked, /DISPLAY, /nest, weight=weight, MASKCUBE=maskcube, masknstacked=masknstacked, mask_stack_output=mask_stack_output

; Create FITS image
  IF KEYWORD_SET(fits) THEN BEGIN
     coord = 'G'                ; can be 'C' for Celestial, 'G' for Galactic and 'E' for Ecliptic
     gnom2fits, img_stack_output, fits, rot=[0., 0.], coord=coord, reso=reso_arcmin
  ENDIF

  PRINT, ' number of sources stacked: ', nstacked

ENDIF
;  STOP


END