;+
; NAME:
;                  AUX_ALONG
; PURPOSE:
;                  Derive age/metallicity profiles for long-slit spectra
; USAGE:
;                  aux_along, <obs>, <cmp>, <target_sn>, <file_out>, /PLOT,
;                        MD=md, KMOMENT=kmoment,
;                        CENTER=center
;
; DESCRIPTION:
;   <obs> shall be a spect structure containing the long-slit spectrum.
;   The program search the peak of light (center of the galaxy),
;   and analyse the data successively from each side of the peak.
;   The data are binned together to reach <target_sn>, but if the 
;   radial range assembled in one bin reach rmax/rmin=1.4 the target_sn
;   is reduced. The profile is stopped when it is not possible to get
;   a bin with sn>5.
;
;   The output of this program is a series of FITS and res files, each
;   of them corresponding to a single point of the profile.
;   The file names are made by attaching 'r'<rad> to the prefix <file_out>,
;   where <rad> is the value of the radius (positive or negative).
;
; ARGUMENTS:
;  <spectrum> -  input 
;                2d spectrum structure (definined with ULY_SPECT_ALLOC) 
;
;  <cmp>      -  input
;                Array of models' components, see ULY_FIT.
;                Prior calls to ULY_STAR, ULY_TGM, ULY_SSP and/or ULY_LINE
;                can define this array and set the guesses and constraints
;                on the free parameters.
;                MODEL_FILE can be given instead.
;
;  <target_sn>-  input
;                Targeted S/N ratio; data are binned together to reach
;                this value
;
;  <file_out> -  input, string
;                output file name, to which suffix with the readious
;                and '.res'/'.fits' are attached
;
; KEYWORDS:
;   /PLOT
;      Plot ulyss fit
;
;    MD 
;      Degree of multiplicative polynomial, to be passed to ulyss (see
;      ulyss.pro documentation)
;
;    KMOMENT 
;      Number of terms in LOSVD, default 2 - velocity, velocity dispersion
;
;    CENTER
;      Index of the center line. By default it is automatically determined
;      as the peak of light.
;
; HISTORY:
;    MK: 2009/09  
;      Creation 
;   PhP: 2010/05/25
;      Do the fit in 2 steps (first clean, after global minimization)
;      Dismiss the bad fits (as much as I could)
;   PhP: 2010/06/23
;      Make the program more general, to be usable for any longslit
;      spectra
;    MK: 2012/10/24 
;      Documentation
;- 
; CATEGORY:
;                  ULY_AUX
  
pro aux_along, obs, cmp, target_sn, file_out, PLOT=plot, MD=md, KMOMENT=kmoment, CENTER=center

  if n_elements(obs) eq 0 then $
     print, 'Usage: AUX_ALONG, obs, cmp, targer_sn, file_out ...'

; We are fixing some parameters ...
     min_snr = 5.
     ad = -1
     quiet = 1                  ; for ulyss

     if n_elements(kmoment) eq 0 then kmoment = 2

     if n_elements(target_sn) eq 0 then message, 'argument <target_sn> is required'

; Find the step along the slit, arcsec/pix
     if n_elements(*obs.hdr) eq 0 then *obs.hdr =  ['']
     step2 = double(string(sxpar(*obs.hdr,'CDELT2',COUNT=count)))
     if count eq 0 then step2 = 1 ; default = 1 arcsec/pix

; Determine the light profile and the position of the center (center)
     mask2d = bytarr(n_elements(*obs.data))
     mask2d[*obs.goodpix] = 1
     dim = size(*obs.data, /DIM)
     mask2d = reform(mask2d, dim)
     goodpix = where(mask2d[*,0] eq 1, cnt)
     if cnt gt 0 then $
        profile = total((*obs.data)[goodpix,*],1,/NaN) $
     else return
                  
     if n_elements(center) eq 0 then begin
        mflux = reform(where(profile eq max(profile))) ; line of max flux
        if not mflux eq 0 then begin 
           n = 4
           m = indgen(2*n+1) - n
           center = total(profile[mflux-n:mflux+n] * m) / (total(profile[mflux-n:mflux+n])-profile[mflux]) + mflux
        endif else $
           center = mflux ;; already centered
     endif
     mflux = nint(center)

     pos = step2 * (findgen(n_elements((*obs.data)[0,*])) - center[0]) ;pos in the slit

; Save the light profile in a FITS file
     fits_open, file_out+'_profile.fits', fcb, /WRITE ; create a new FITS file
     h = uly_spect_get(obs, /HDR)
     if n_elements(h) le 1 then h = ['']
     sxaddpar, h, 'CRVAL1', 0.
     sxaddpar, h, 'CRPIX1', center[0], 'position of the center'
     sxaddpar, h, 'TITLE', 'light profile along the slit'
     sxaddpar, h, 'CDELT1', step2
     fits_write, fcb, profile, h
     fits_close, fcb

     obs0 = uly_spect_alloc(SPECTRUM=obs)

; save the guesses
     if ptr_valid((*cmp[0].para)[0].guess ) then begin
        guess0 = *(*cmp[0].para)[0].guess
        guess1 = *(*cmp[0].para)[1].guess
     endif
;; re-init the guesses
;*(*cmp[0].para)[0].guess = alog([1000d,4000d]) ;  AG=[1000.,4000]
;*(*cmp[0].para)[1].guess = [-1d,-0.35d,0d] ; zg=[-1.,-0.35,0]

;; for the right side
     t_sn = target_sn
     i = mflux
     mm = 0
     nscan = n_elements((*obs0.data)[0,*])
     while i ge mflux and i+mm+1 lt n_elements(pos)-1 do begin
        mm = 0
        ma = 0
        data = (*obs.data)[*,i]
        err2 = (*obs.err)[*,i]^2
        mask = mask2d[*,i]
        med_snr = median(data/sqrt(err2))
        if finite(med_snr) eq 0 then begin
           med_snr = 0
           data = 0
           err2 = 0
        endif
        if  med_snr lt t_sn then begin      
           while  med_snr lt t_sn and i+mm lt nscan - 1 do begin
              mm += 1
              data_n = data + (*obs.data)[*,i+mm]
              err2_n = err2 + (*obs.err)[*,i+mm]^2
              mask_n = mask + mask2d[*,i+mm]
              snr = median(data_n/sqrt(err2_n))
              
              if finite(snr) eq 0 then snr = -1
              if snr ge med_snr then begin ; keep only good data
                 data = data_n
                 err2 = err2_n
                 mask = mask_n
                 med_snr = snr
                 ma = mm
              endif
              if mm gt 5 and pos[i+ma]/pos[i] gt 1.4 and t_sn gt min_snr then begin
                 t_sn = max([med_snr, min_snr])
                 print, 'reduce the targeted SNR to', t_sn
              endif
           endwhile
        endif
        if med_snr lt min_snr then break ; too low to analyse

        posi = (pos[i] + pos[i+ma]) / 2d
        print, 'pos:'+strtrim(pos[i+ma/2],2) + ' snr:' + strtrim(med_snr,2)

        *obs0.data = data
        *obs0.err = sqrt(err2)
        *obs0.goodpix = where(mask eq 1b+mm, goodcnt) ;; 1 good!
        if goodcnt le 1 then break  
        neg = where (*obs0.data/*obs0.err le 0.01*med_snr, cnt)
        if cnt gt 0 then begin
           mask = uly_spect_get(obs0, /MASK)
           mask[neg] = 0
           plus = where (mask gt 0, cnt)
           if cnt gt 0 then *obs0.goodpix = plus
        endif

        spi = where (*obs0.data ge 10*median(*obs0.data), cnt)
        if cnt gt 0 then begin
           print, 'spikes:', cnt, median(*obs0.data)
           mask = uly_spect_get(obs0, /MASK)
           plus = where (mask gt 0, cntp)
           mask[spi] = 0
           plus = where (mask gt 0, cnt)
           if cnt gt 0 then *obs0.goodpix = plus
        endif

        weight = 1d / err2
        large_weight = where(weight gt 100*mean(weight), cnt, COMP=small_weight)
        if cnt gt 0 then begin
           meanerr = sqrt(1d/mean(weight[small_weight]))
           (*obs0.err)[large_weight] = meanerr / 4
        endif

        file = file_out+'_r'+strtrim(posi,2)
        ulyss, obs0, cmp, FILE_OUT=file, AD=ad, /CLEAN, MD=md, QUIET=quiet, KMOMENT=kmoment
        undefine, str
        if file_test(file+'.res') then file_delete,file+'.res'
        if file_test(file+'.fits') then ulyss, file+'.fits', cmp, FILE_OUT=file, PLOT=plot, AD=ad, MD=md, QUIET=quiet, SOLUT=str, KMOMENT=kmoment
;    ulyss, obs0, cmp, FILE_OUT=file, PLOT=plot, AD=ad, /CLEAN, MD=md, QUIET=quiet, SOLUT=str

; search if any of the params reached its limit
        if size(str,/TYPE) eq 8 then begin
           ok = 1
           if ptr_valid((*cmp[0].para)[0].guess ) then begin ;; only if we fit a cmp with parameters
              if (*str[0].cmp[0].para)[0].value eq (*str[0].cmp[0].para)[0].limits[0] then begin
                 print, 'parameter 0 reached its lower limits'
                 ok = 0
              endif
              if (*str[0].cmp[0].para)[0].value eq (*str[0].cmp[0].para)[0].limits[1] then begin
                 print, 'parameter 0 reached its upper limits'
;            ok = 0   Do not reject this case !!!
              endif
              if (*str[0].cmp[0].para)[1].value eq (*str[0].cmp[0].para)[1].limits[0] or $
                 (*str[0].cmp[0].para)[1].value eq (*str[0].cmp[0].para)[1].limits[1] then begin
                 print, 'parameter 1 reached one of the limits'
                 ok = 0
              endif
              if (*str[0].cmp[0].para)[0].fixed eq 0 then begin
                 if (*str[0].cmp[0].para)[0].error eq 0 then begin
                    print, 'parameter 0 has a null error' 
;            ok = 0  Do not reject !!!
                 endif else if (*str[0].cmp[0].para)[0].error le 0.001 then begin
                    print, 'parameter 0 has a small error' 
                 endif
              endif
              if (*str[0].cmp[0].para)[1].fixed eq 0 then begin
                 if (*str[0].cmp[0].para)[1].error eq 0 then begin
                    print, 'parameter 1 has a null error' 
                    ok = 0
                 endif else if (*str[0].cmp[0].para)[1].error le 0.001 then begin
                    print, 'parameter 1 has a small error' 
                 endif
              endif 
           endif
           if str[0].e_losvd[1] gt 5.*str[0].losvd[1] then begin
              print, 'velocity dispersion has a too large error' 
              ok = 0
           endif
           if ok eq 0 then begin
              file_delete, file + '.res'
              file = file_out + '_REJECT' +strtrim(posi,2)
              uly_solut_twrite, str, file
           endif
        endif

        i += mm + 1

;  for next scan use the present age/met as guess
;    str = uly_solut_tread(file)
        if size(str,/TYPE) eq 8 then begin
           if ptr_valid((*cmp[0].para)[0].guess ) then begin
              age = (*str[0].cmp[0].para)[0].value
              met = (*str[0].cmp[0].para)[1].value
              if age lt alog(8000d) and age ge alog(100d) then $
                 *(*cmp[0].para)[0].guess = age - [0, alog(1.5)]
              if (*cmp[0].para)[1].fixed eq 0 and met lt 0.3 and met ge -1.8 then $
                 *(*cmp[0].para)[1].guess = met + [0d, -0.2d]
              heap_free, str
           endif
        endif

     endwhile

; reinit the guess
     if ptr_valid((*cmp[0].para)[0].guess ) then begin
        *(*cmp[0].para)[0].guess = guess0
        *(*cmp[0].para)[1].guess = guess1
     endif
; for the left side
     t_sn = target_sn
     i = mflux-1
     mm = 0
     ma = 0
     while i le mflux and i-mm-1 gt 0 do begin
        mm = 0
        data = (*obs.data)[*,i]
        err2 = (*obs.err)[*,i]^2
        mask = mask2d[*,i]
        med_snr = median(data/sqrt(err2))
        if finite(med_snr) eq 0 then begin
           med_snr = 0
           data = 0
           err2 = 0
        endif
        if  med_snr lt t_sn then begin      
           while  med_snr lt t_sn and i-mm gt 0 do begin
              mm += 1
              data_n = data + (*obs.data)[*,i-mm]
              err2_n = err2 + (*obs.err)[*,i-mm]^2
              mask_n = mask + mask2d[*,i-mm]
              snr = median(data_n/sqrt(err2_n))
              if finite(snr) eq 0 then snr = -1
              if snr ge med_snr then begin ; keep only good data
                 data = data_n
                 err2 = err2_n
                 mask = mask_n
                 med_snr = snr
                 ma = mm
              endif
              if mm gt 5 and pos[i-ma]/pos[i] gt 1.4 and t_sn gt min_snr then begin
                 t_sn = max([med_snr, min_snr])
                 print, 'reduce the targeted SNR to', t_sn
               endif
           endwhile
        endif
        *obs0.data = data
        *obs0.err = sqrt(err2)
        *obs0.goodpix = where(mask eq 1b+mm, goodcnt) ;; 1 good!

        if goodcnt le 1 then break

        if med_snr lt min_snr then break ; too low to analyse

        posi = (pos[i] + pos[i-ma]) / 2d
        print, 'pos:'+strtrim(posi,2) + ' snr:' + strtrim(med_snr,2)
        neg = where (*obs0.data/*obs0.err le 0.01*med_snr, cnt)

        if cnt gt 0 then begin
           mask = uly_spect_get(obs0, /MASK)
           mask[neg] = 0
           plus = where (mask gt 0, cnt)
           if cnt gt 0 then *obs0.goodpix = plus
        endif

        spi = where (*obs0.data ge 10*median(*obs0.data), cnt)
        if cnt gt 0 then begin
           print, 'spikes:', cnt, median(*obs0.data)
           mask = uly_spect_get(obs0, /MASK)
           mask[spi] = 0
           plus = where (mask gt 0, cnt)
           if cnt gt 0 then *obs0.goodpix = plus
        endif

        weight = 1d / err2
        large_weight = where(weight gt 100*mean(weight), cnt, COMP=small_weight)
        if cnt gt 0 then begin
           meanerr = sqrt(1d/mean(weight[small_weight]))
           (*obs0.err)[large_weight] = meanerr / 4
        endif

        file = file_out+'_r'+strtrim(posi,2)

        ulyss, obs0, cmp, FILE_OUT=file, AD=ad, /CLEAN, MD=md, QUIET=quiet, KMOMENT=kmoment
        undefine, str
        if file_test(file+'.res') then file_delete,file+'.res'
        if file_test(file+'.fits') then ulyss, file+'.fits', cmp, FILE_OUT=file, PLOT=plot, AD=ad, MD=md, QUIET=quiet, SOLUT=str, KMOMENT=kmoment

; search if any of the params reached its limit
        if size(str,/TYPE) eq 8 then begin
           ok = 1
           if (*str[0].cmp[0].para)[0].value eq (*str[0].cmp[0].para)[0].limits[0] then begin
              print, 'parameter 0 reached its lower limits'
              ok = 0
           endif
           if (*str[0].cmp[0].para)[0].value eq (*str[0].cmp[0].para)[0].limits[1] then begin
              print, 'parameter 0 reached its upper limits'
;            ok = 0   Do not reject this case !!!
           endif
           if (*str[0].cmp[0].para)[1].value eq (*str[0].cmp[0].para)[1].limits[0] or $
              (*str[0].cmp[0].para)[1].value eq (*str[0].cmp[0].para)[1].limits[1] then begin
              print, 'parameter 1 reached one of the limits'
              ok = 0
           endif
           if (*str[0].cmp[0].para)[0].error eq 0 then begin
              print, 'parameter 0 has a null error' 
;            ok = 0  Do not reject !!!
           endif else if (*str[0].cmp[0].para)[0].error le 0.001 then begin
              print, 'parameter 0 has a small error' 
           endif
           if (*str[0].cmp[0].para)[1].error eq 0 then begin
              print, 'parameter 1 has a null error' 
              ok = 0
           endif else if (*str[0].cmp[0].para)[1].error le 0.001 then begin
              print, 'parameter 1 has a small error' 
           endif
           if str[0].e_losvd[1] gt 5.*str[0].losvd[1] then begin
              print, 'velocity dispersion (',str[0].losvd[1],')has a too large error (',str[0].e_losvd[1],')'
              ok = 0
           endif
           if ok eq 0 then begin
              file_delete, file + '.res'
              file = file_out + '_REJECT' +strtrim(posi,2)
              uly_solut_twrite, str, file
           endif
        endif

        i -= (mm + 1)

;  for next scan use the present age/met as guess
;    str = uly_solut_tread(file)
        if size(str,/TYPE) eq 8 then begin
           if ptr_valid((*cmp[0].para)[0].guess ) then begin
              age = (*str[0].cmp[0].para)[0].value
              met = (*str[0].cmp[0].para)[1].value
              if age lt alog(8000d) and age ge alog(100d) then *(*cmp[0].para)[0].guess = age - [0, alog(1.5)]
              if (*cmp[0].para)[1].fixed eq 0 and met lt 0.3 and met gt -1.8 then *(*cmp[0].para)[1].guess = met + [0d, -0.2d]
              heap_free, str
           endif
        endif

     endwhile

; reinit the guess
     if ptr_valid((*cmp[0].para)[0].guess ) then begin
        *(*cmp[0].para)[0].guess = guess0
        *(*cmp[0].para)[1].guess = guess1
     endif
     uly_spect_free, obs0

  end
;------------- end -----------------