;; input a file with RA, DEC fields.  
;; this script assumes the following fields already exist in the input structure, to be filled in with this code:
;; AIRMASS, EB_MINUS_V, SKY_FLUX, PSF_FWHM
;; imagingdropt: 0 is just straight DR8 imaging., 1 is go by chunk geometry.

;; call as fill_BOSS_imaging,
;; fill_BOSS_imaging,"../data/trimmed-collate-cmasspy-dr12v0",0,0,/verboseopt
;; tagged that output with _dr8imaging, now run with imagingdropt = 1
;; fill_BOSS_imaging,"../data/trimmed-collate-cmasspy-dr12v0",1,0,/verboseopt


pro fill_BOSS_imaging,infilebase,imagingdropt,DorRopt,verboseopt=verboseopt

  ;; first just use DR8 imaging for this.
  ;; extra fields will be added here, so I can take them out of mksample I guess.  Let's do that after testing.
  ;; for data we'll use the ichunk field from the collate file.
  ;; for the randoms we'll use chunk geometry files.  
  ;; some boundaries will get screwed up, I guess.

  ;; imagingdropt = 0; this is for final imaging only.

  ;; setup chunk geometry stuff.  copying from gather_photometry_OLD.pro
  start_chunk = 2
  end_chunk   = 15      ; Beyond this photometry is "frozen"; resolve is frozen past 11
  nchunk      = end_chunk-start_chunk + 1
  ; Photo_tag stores the version of the photometry to be used, and
  ; whether the file is already opened by the code.
  ; Column 0: chunk number.
  ; Column 1: photo version.
  ; Column 2: whether photo file already opened.
  photo_tag_all      = lonarr(nchunk,3)
  photo_tag_all[*,0] = indgen(nchunk) + start_chunk ; chunk number 
  for ichunk=0,nchunk-1L do begin
    i = photo_tag_all[ichunk,0]
    if   (i EQ 1 )                  then  photo_tag_all[ichunk,1] = 0 ; comm
    if   (i EQ 2 )                  then  photo_tag_all[ichunk,1] = 1 ; comm2
    if ( (i GE  3) AND ( i LE  4) ) then  photo_tag_all[ichunk,1] = 2 ; main002
    if ( (i GE  5) AND ( i LE  6) ) then  photo_tag_all[ichunk,1] = 5 ; main005
    if ( (i GE  7) AND ( i LE 11) ) then  photo_tag_all[ichunk,1] = 7 ; main007
    if ( (i GE 12) AND ( i LE 14) ) then  photo_tag_all[ichunk,1] = 8 ; main008
    if   (i GE 15)                  then  photo_tag_all[ichunk,1] = 9 ; main009
    ; for galaxies this is frozen from chunk15 and onwards
  endfor
  photo_tag_all[*,2] = -1

  dir_bosstilelist= filepath(root=getenv("BOSSTILELIST_DIR"),"/")
  ;;; end setup chunk geometry stuff.

  ; Set up some useful paths for the code below.
  data_dir  = filepath(root=getenv("BOSSLSS_DIR"),"/data/")
  
  inoutfile = data_dir + infilebase+'.fits'
  print, 'file name is ',inoutfile
  dd =  hogg_mrdfits(inoutfile,1,hdr,nchunk=50)

  print,'MAKE SURE THIS IS CORRECTLY SET TO BOSS!!'
  ;; riemann setting.
  ;setenv,'PHOTO_RESOLVE=/clusterfs/riemann/raid006/bosswork/boss/resolve/2010-05-23'
  ;; generalize to any system where PHOTO_RESOLVE is set.
  a = strsplit(getenv('PHOTO_RESOLVE'), "resolve", /extract, /regex)
  resolvebase = a[0] + 'resolve/'
  newp = resolvebase + '2010-05-23'
  setenv,'PHOTO_RESOLVE='+newp
  print,getenv('PHOTO_RESOLVE')

  ;; first use final DR8 imaging and fill in everything.
  find_primary_field,dd.ra,dd.dec,fout=fout,ilist=ilist,select_tags=['run','rerun','camcol','field','airmass','skyflux','psf_fwhm','image_depth'],/unified

  dd.airmass = fout.airmass
  dd.psf_fwhm = fout.psf_fwhm
  dd.skyflux = fout.skyflux
  dd.image_depth = fout.image_depth

  ;; hard-coding resolve versions from here:
  ;; https://trac.sdss3.org/wiki/BOSS/target
  if imagingdropt eq 1 then begin
    ;; go back and fill in different values when earlier photo_resolve reductions were used for targeting, using chunk geometry files.
    ;;for ichunk=2L,11L do begin
    ;; boss-gal-clustering email from Jeremy Tinker on June 4 2014, lower chunk number has priority.
    ;; I ran the data through this procedure in chunktest.pro and chunktestS.pro (see chunktestout)
    ;; I found excellent (basically perfect) agreement with the ichunk field for the galaxies.
    for ichunk=11L,2L,-1L do begin
      print,'working on chunk',ichunk 
      if ichunk ge 1 and ichunk le 4 then begin
        newp = resolvebase + '2009-06-14'
        setenv,'PHOTO_RESOLVE='+newp
        ; path on riemann.
        ;setenv,'PHOTO_RESOLVE=/clusterfs/riemann/raid006/bosswork/boss/resolve/2009-06-14'
      endif
      if ichunk ge 5 and ichunk le 11 then begin
        newp = resolvebase + '2009-11-16'
        setenv,'PHOTO_RESOLVE='+newp
        ; path on riemann.
        ;setenv,'PHOTO_RESOLVE=/clusterfs/riemann/raid006/bosswork/boss/resolve/2009-11-16'
      endif
      print,getenv('PHOTO_RESOLVE'),'used for ichunk = ',ichunk

      if DorRopt EQ 0 then begin ;; data.  just use ichunk field in the collate file to select them.
        xxchunk = where(dd.ichunk eq ichunk,cntxx)      
      endif else begin  ;; see difference between ichunk and geometric output in chunktest.21579.riemann.out
        chunkname = 'boss'+strtrim(string(ichunk),2)
        polyname  = dir_bosstilelist+'/outputs/'+chunkname+$
                    '/geometry-'+chunkname+'.ply'
        read_mangle_polygons,polyname,chunk_geom
        ; Find the objects that are in the current chunk
        xxchunk = where(is_in_window(ra=dd.ra,dec=dd.dec,chunk_geom)$
           GT 0,cntxx)
      endelse
      if cntxx GT 0 then begin
        print,'reassigning this many elements to ichunk ',ichunk,n_elements(xxchunk)
        ;; unified gives an error:
        ;; % FIND_PRIMARY_FIELD1: Unified field list missing from $PHOTO_RESOLVE.  Just run the slow way instead.
        ;; works now??
        find_primary_field,dd[xxchunk].ra,dd[xxchunk].dec,fout=fouttmp,ilist=ilisttmp,select_tags=['run','rerun','camcol','field','airmass','skyflux','psf_fwhm','image_depth'],/unified
        ;;find_primary_field,dd[xxchunk].ra,dd[xxchunk].dec,fout=fouttmp,ilist=ilisttmp,select_tags=['run','rerun','camcol','field','airmass','skyflux','psf_fwhm','image_depth']
        dd[xxchunk].airmass = fouttmp.airmass
        dd[xxchunk].psf_fwhm = fouttmp.psf_fwhm
        dd[xxchunk].skyflux = fouttmp.skyflux
        dd[xxchunk].image_depth = fouttmp.image_depth
        ;; fout[xxchunk] = fouttmp
        if n_elements(fouttmp) ne n_elements(fout[xxchunk]) then begin
          print,'no good!',n_elements(fouttmp),n_elements(fout[xxchunk])
          stop
        endif
        ;struct_assign,fouttmp,fout[xxchunk] ;; did this really work? No, error.
        ;struct_assign,ilisttmp,ilist[xxchunk] ;; did this really work?
        fout[xxchunk].run = fouttmp.run
        fout[xxchunk].rerun = fouttmp.rerun
        fout[xxchunk].camcol = fouttmp.camcol
        fout[xxchunk].field = fouttmp.field
        fout[xxchunk].airmass = fouttmp.airmass
        fout[xxchunk].skyflux = fouttmp.skyflux
        fout[xxchunk].psf_fwhm = fouttmp.psf_fwhm
        fout[xxchunk].image_depth = fouttmp.image_depth
        ilist[xxchunk] = ilisttmp

      endif ;; have items in the chunk.
    endfor ;; ichunk loop.

  endif
  if keyword_set(verboseopt) then begin
    print,'verbose is set!'
    outfile1 = infilebase+'_fout.fits'
    outfile2 = infilebase+'_ilist.fits'
    ;; will compare with original in separate python script.
    mwrfits,fout,outfile1,/create
    mwrfits,ilist,outfile2,/create
  endif

  ;; finally, fill in eb_minus_v
  ;; jeremy used this function, not euler.  which is correct?  i noticed ~4 armin differences (max) with euler,
  ;; though vast majority were in 2'' agreement.
  ;; glactc, dd.ra, dd.dec, 2000., gl, gb, 1, /deg
  euler,dd.ra,dd.dec,gl,gb,1
  dd.eb_minus_v = dust_getval(gl,gb, /interp, /noloop)

  mwrfits,dd,inoutfile,/create
  

;  stop
end