function g_r2age,g_r
    coeff=[-1.1157153, 3.2469886, -5.2583179, 5.2132806]
    age=10^poly(g_r,coeff)
    return,(0.1 > age < 12.0)
end

function gr_from_sdsshdr,h,ebv=ebv
    ebv=sxpar(h,'SFD_EBV')
    mag_g=sxpar(h,'MAG_G')-ebv*3.793
    mag_r=sxpar(h,'MAG_R')-ebv*2.751
    z=sxpar(h,'Z')
    kcgc=[0.00000,0.00000,0.00000,0.00000,-0.965477,2.26834,4.25947,-3.39657,15.7355,-15.5160,20.3434,4.50132,-81.8360,-52.1578,-43.2650,275.658,139.375,-321.553]
    kcg=poly2d(z,mag_g-mag_r,kcgc,/irr,deg1=5,deg2=3)
    kcrc=[0.00000,0.00000,0.00000,0.00000,-0.366693,2.67325,-2.93288,1.34941,1.86297,15.3022,0.568358,-3.81143,-68.5488,-58.5031,17.6783,263.468,31.5546,-258.72]
    kcr=poly2d(z,mag_g-mag_r,kcrc,/irr,deg1=5,deg2=3)
    return,[mag_g-kcg,mag_r-kcr]
end

pro process_sdss, inptable=inptable, outpath=outpath, datapath=datapath,$
    plot=plot, moments=moments, mdegree=mdegree, degree=degree,$
    lammin=lammin, lammax=lammax, tmpdir=tmpdir, nbursts=nbursts,$
    emexcl=emexcl, start=start, fakesig=fakesig, corrfakesig=corrfakesig, $
    iterate=iterate, maxchi2=maxchi2, path_ssp=path_ssp,suffix=suffix,$
    bspline=bspline,miles=miles,expsfh=expsfh,$
    inpmetidx=inpmetidx,inpageidx=inpageidx,fixpar=fixpar,imin=imin,imax=imax,$
    phweight=phweight,sed_data=sed_data,sed_errors=sed_errors,nophoterror=nophoterror,$
    sed_path=sed_path,sed_suffix=sed_suffix,allphotbands=allphotbands,goodbands=goodbands,sed_chi2fld=sed_chi2fld,$
    nlosvd=nlosvd,bias=bias,broademl=broademl

age_universe=13750.0
maxage=19000.0
if(n_elements(phweight) ne 1) then phweight=0d

if(n_elements(maxchi2) ne 1) then maxchi2=0.8
if(n_elements(nbursts) ne 1) then nbursts=1
if(n_elements(nlosvd) ne 1) then nlosvd=1
if(n_elements(fixpar) ne nlosvd*2+nbursts*2) then fixpar=intarr(nlosvd*2+nbursts*2)
if(n_elements(tmpdir) ne 1) then begin
    tmpdir=string(randomu(seed,/long),format='(i20.20)')+'sdss'
    spawn,'mkdir -p '+tmpdir
    deltmp=1
endif else deltmp=0
if(n_elements(mdegree) ne 1) then mdegree=5
if(n_elements(moments) ne 1) then moments=2
if(n_elements(inptable) ne 1) then inptable='tabsdss.txt'
if(n_elements(outpath) ne 1) then outpath='SDSS/results/'
if(n_elements(datapath) ne 1) then $
     datapath='../processed/SDSS/data/'
if(n_elements(lammin) ne 1) then lammin=3910.0
if(n_elements(lammax) ne 1) then lammax=6790.0

if(keyword_set(miles)) then begin
    if(n_elements(path_ssp) eq 0) then $
    path_ssp='/usr/local/rsi/idl71/external/mod/templates/template/pegase.hr/MILES/' ;;; PEGASE-MILES
stop
;;;;    path_ssp='../template/Vazdekis/KU/output/' ;;; Original Vazdekis models
;;;;    path_ssp='../template/pegase.hr/MILES/' ;;; PEGASE-MILES
    suffix='_KU.fits'
;;;    inpageidx=range(0,49)
    inpageidx=range(0,67)
    inpmetidx=range(1,7)
endif

tab=read_asc(inptable,/str)
s_tab=size(tab)
n_gal=s_tab[2]
if(s_tab[0] eq 1) then n_gal=1

if(n_elements(imin) ne 1) then imin=0L
if(n_elements(imax) ne 1) then imax=n_gal-1L

for i=imin,imax do begin
  spData = readfits(tab[0,i], hpri)
  dirsuff=(s_tab[1] eq 2)? tab[1,i]+'/' : '/'
  parse_spechdr,hpri,wl=wlvac,velScale=velScaleSDSS
  wlatm=wlvac/(1.0d + 2.735182d-4 + 131.4182/wlvac^2 + 2.76249d+8/wlvac^4)
  linData = mrdfits(tab[0,i], 2, hlin)
  resData = mrdfits(tab[0,i], 6, hres)
  z_sdss = sxpar(hpri,'Z')
  sig_sdss = sxpar(hpri,'VEL_DIS')
  if(sig_sdss lt 10.0) then sig_sdss=100.0
  print,'Processing ',i,'/',imax,' (',i-imin+1L,'/',imax-imin+1L,')',format='(a,i7,a,i7,a,i7,a,i7,a)'
  print,'VELSCALE,V,SIG=',velscalesdss,z_sdss*299792.5d,sig_sdss

  oiiiline = linData[30]
  haline = linData[36]
  niiline = linData[37]
  emlines = (((oiiiline.ew gt 0.) or (haline.ew gt 0.) or (niiline.ew gt 0.)) $
         or keyword_set(emexcl))? 1 : 0;
  if(n_elements(emexcl) eq 1) then if(emexcl eq 0) then emlines=0
  basename = strcompress(string(sxpar(hpri,'MJD'),format='(i5.5)')+$
                     '-'+string(sxpar(hpri,'PLATEID'),format='(i4.4)')+$
                     '-'+string(sxpar(hpri,'FIBERID'),format='(i3.3)'),$
                     /remove_all)
  tmpname = 'spSpecResmp-'+basename+'.fits'
  etmpname = 'err_'+tmpname
  sxaddpar,hpri,'CTYPE1','WAVE-LOG'
  sxaddpar,hpri,'CDELT1',sxpar(hpri,'CD1_1')
  mederr=median(spData[*,2])
  bdata=where((spData[*,2] le 0.001) or (spData[*,2] gt mederr*4e+4), cbdata)
  if(cbdata gt 0) then spData[bdata,*]=!values.f_nan
  writefits,tmpdir+'/'+tmpname,spData[*,0],hpri
  writefits,tmpdir+'/'+etmpname,spData[*,2],hpri

  writefits,tmpdir+'/'+'velmap.fts',double(((wlvac-wlatm)/wlatm)*299792.458)

  if(keyword_set(miles)) then begin
    gkrn = psf_gaussian(Npixel=31,fwhm=[2.355*80.0/velScaleSDSS],/norm,ndim=1)
    ;writefits,tmpdir+'/'+tmpname,convol(spData[*,0],gkrn),hpri
    ;writefits,tmpdir+'/'+etmpname,convol(spData[*,2],gkrn),hpri
    c_disp=[53.193741,-1.0550717,88.539196,-222.85372,$
            -806.58457,995.97896,3327.8268,-2748.5269,$
            -6917.4259,4944.2669,7173.6642,-5076.2970,-2950.6301,2181.4498]
    lam_cnt=5193.6121
    lam_norm=2000.0
    disp_corr=poly((wlvac-lam_cnt)/lam_norm,c_disp)
    b_wl=where(wlvac gt 7200,cb_wl)
    if(cb_wl gt 0) then disp_corr[b_wl]=disp_corr[b_wl[0]]
    disp_vec=sqrt(0*80d^2+(double(resData.DISPERSION)*velScaleSDSS)^2-disp_corr^2)
    bdv=where(finite(disp_vec) ne 1 or disp_vec le 10.0,cbdv)
    if(cbdv gt 0) then disp_vec[bdv]=10.0
    writefits,tmpdir+'/'+'dismap.fts',disp_vec
  endif else begin
    disp_vec=double(resData.DISPERSION)*velScaleSDSS
    bdv=where(finite(disp_vec) ne 1 or disp_vec le 10.0,cbdv,compl=gdv)
    if(cbdv gt 0) then disp_vec[bdv]=interpol(disp_vec[gdv],wlvac[gdv],wlvac[bdv])
    bdv=where(finite(disp_vec) ne 1 or disp_vec le 10.0,cbdv,compl=gdv)
    if(cbdv gt 0) then disp_vec[bdv]=10.0
    writefits,tmpdir+'/'+'dismap.fts',disp_vec
  endelse

  if(keyword_set(expsfh)) then begin
      maxage=50000.0
      ;path_ssp='../template/pegase.hr/grid_expSFH_kroupa/HR/'
      path_ssp='/usr/local/rsi/idl71/external/mod/templates/template/pegase.hr/grid_expSFH_kroupa/HR/'
      gal_age=galage((z_sdss > 0.0),1000.0,h0=72.0,/silent)/1e+6 ;;; z_form=1000.0
      age_min = (2 > (fix(gal_age/1000.0)) < 13)*1000
      age_max = age_min+1000
      w_min = 1.0 - ((0.0 > ((gal_age - age_min)/1000.0) < 1.0))
      tmpgridweights=[w_min, 1.0-w_min]
      
      suffix=['_expall_age'+string(age_min,format='(i5.5)')+'Myr_Kroupa_0.1_120_elodie31.fits',$
              '_expall_age'+string(age_max,format='(i5.5)')+'Myr_Kroupa_0.1_120_elodie31.fits']

      ;print,gal_age,suffix,tmpgridweights
      inpageidx=range(0,11)
      inpmetidx=range(1,9)
  endif

  writefits,tmpdir+'/'+'h3map.fts',wlvac*0d
  writefits,tmpdir+'/'+'h4map.fts',wlvac*0d
  writefits,tmpdir+'/'+'wl.fts',wlvac
;;  if((n_elements(start) ne nbursts*2+2) or (not(arg_present(start)))) then begin
  if((n_elements(start) ne nbursts*2+2*nlosvd)) then begin
      gr_sdss=gr_from_sdsshdr(hpri)
      age_st=g_r2age(gr_sdss[0]-gr_sdss[1])*1000.0
      if(keyword_set(expsfh)) then age_st=(20.0 > (13000.0-age_st) < 8000.0)
      start=[z_sdss*299792.458,sig_sdss,age_st,-0.2]
      for hh=1,nbursts-1 do start=[start,2000.0/double(hh),0.0]
      for hh=1,nlosvd-1 do start=[z_sdss*299792.458,sig_sdss,start]
  endif else start[0:1]=[z_sdss*299792.458,sig_sdss]
  z=z_sdss
  ;exclreg=(emlines eq 1)? [[4861.*(1.+z),3.0],[4959.*(1.+z),3.0],[5007.*(1.+z),3.0],[5577.0,4.0],[5645.0,4.0]] : [[5577.0,4.0],[5645.0,4.0]]
  eml_exclreg=[[3727.0*(1.+z),10.0],[3698.0*(1.+z),7.0],[3703.0*(1.+z),7.0],[3711.0*(1.+z),7.0],[3750.0*(1.+z),7.0],[3770.0*(1.+z),7.0],[3798.0*(1.+z),7.0],[3836.0*(1.+z),7.0],[3868.0*(1.+z),7.0],[3889.0*(1.+z),7.0],[3970.0*(1.+z),6.0],[4103.0*(1.+z),7.0],[4341.0*(1.+z),9.0],[4364.0*(1.+z),8.0],[4472.0*(1.+z),8.0],[4687.0*(1.+z),9.0],[4862.7*(1.+z),12.0],[4960.3*(1.+z),12.0],[5008.2*(1.+z),13.0],[5200.*(1.+z),7.],[5877.*(1.+z),10.],[6301.*(1.+z),7.],[6364.*(1.+z),7.],[6550.*(1.+z),12.],[6564.*(1.+z),14.],[6585.*(1.+z),14.],[6678.*(1.+z),8],[6719.*(1.+z),10.],[6732.*(1.+z),10.],[5894.5*(1.+z),14.0],[5578.0,8.0],[5646.0,8.0],[5894.5,14.0],[6300.0,8.0],[6363.0,8.0]]
  exclreg=(emlines eq 1)? eml_exclreg : [[5578.0,8.0],[5646.0,8.0],[5894.5,14.0],[5894.5*(1.0+z),14.0],[6300.0,8.0],[6363.0,8.0]]
  if(emlines eq 1) then begin
    message,/inf,'Masking emission lines'
    n_em_l=n_elements(eml_exclreg[0,*])
    if(n_elements(broademl) eq 1) then exclreg[1,0:n_em_l-7]=exclreg[1,0:n_em_l-7]*broademl
  endif
  f=tmpdir+'/'+tmpname
  ef=tmpdir+'/'+etmpname
  if(phweight gt 0) then begin
    ph_f=tmpdir+'/mags.fits'
    ph_ef=tmpdir+'/errmags.fits'
    writefits,ph_f,sed_data[*,i]
    writefits,ph_ef,sed_errors[*,i]
  endif
  message,/info,'Processing: '+tab[0,i]
  h_orig=headfits(f)
  gp=[range(5,3800)]
  print,'START=',start

  iter=1
  for iternum=0,1 do begin
    if iter eq 0 then continue
    ppxf_1d3d,galfile=f,errfile=ef,$
        velf=velf,disf=disf,h3f=h3f,h4f=h4f,chi2f=chi2f,$
        agef=agef,metf=metf,ampf=ampf,start=start,plot=plot,$
        mom=moments,goodPixels=gp,mdeg=mdegree,degree=degree,$
        lmin=((lammin*(1.+z)) > 3796.0),lmax=(lammax*(1.+z) < 9150.0),calibpath=tmpdir,$
        nbursts=nbursts,$
        exclreg=exclreg,fixpar=fixpar,$
        galbin=galbin,errbin=errbin,bestfitarr=bestfit,wl_new=wl_new,$
        lammin=((lammin*(1.+z)) > 3796.0),lammax=(lammax*(1.+z) < 9150.0),velscale=velScaleSDSS,$
        velbin=velbin,disbin=disbin,fakesig=fakesig,corrfakesig=corrfakesig,$
        path_ssp=path_ssp,suffix=suffix,$
        bspline=bspline,inpageidx=inpageidx,inpmetidx=inpmetidx,tmpgridweights=tmpgridweights,$
        phweight=phweight,photgalfile=ph_f,photerrfile=ph_ef,nophoterror=nophoterror,$
        sed_path=sed_path,sed_suffix=sed_suffix,allphotbands=allphotbands,goodbands=goodbands,sed_chi2fld=sed_chi2fld,nlosvd=nlosvd,$
        bias=bias,flaggoodpix=flagGoodPix,allcompbestfitarr=allcompbestfitarr
  
    iter=0
    if(((chi2f[0] ge maxchi2) or (agef[0] ge maxage)) and keyword_set(iterate)) then begin
        iter=1
        if(chi2f[0] ge maxchi2) then begin
            emlines=1
            exclreg=eml_exclreg
            exclreg[1,*]=exclreg[1,*]*1.3
        endif
        if(nbursts gt 1) then $
            start[findgen(nbursts)*2+2*nlosvd]=((start[findgen(nbursts)*2+2*nlosvd]/3.5) > 250.0)
        print,'Refitting: emlines=',emlines,' START=',start
    endif
  endfor

  if(nlosvd eq 1 and nbursts eq 1) then allcompbestfitarr=!values.f_nan
  new_gal=dblarr(n_elements(galbin),3)
  new_gal[*,0]=galbin
  new_gal[*,1]=errbin
  new_gal[*,2]=bestfit
  
  h_new=h_orig
  sxaddpar,h_new,'CTYPE1','WAVE-ALOG',' logarithmical binning'
  sxaddpar,h_new,'CRVAL1',alog10(wl_new[0]),' log10(wl0)'
  sxaddpar,h_new,'CRPIX1',1,' reference pixel'
  sxaddpar,h_new,'CDELT1',velScalesdss / 299792.458d/alog(10d),' logarithmic step'
  sxaddpar,h_new,'CD1_1',velScalesdss/299792.458d/alog(10d) ,' logarithmic step'

  out_hdr_res=["XTENSION= 'BINTABLE'    / this is the binary table extension"]
  sxaddpar,out_hdr_res,'EXTNAME','SPECTRUM','name of this binary table extension'

  mod_type = (keyword_set(expsfh))? 'EXPSFH' : 'SSP'
  sxaddpar,out_hdr_res,'MOD_TYPE',mod_type
  if(keyword_set(expsfh)) then begin
    sxaddpar,out_hdr_res,'AGESTART',gal_age
  endif

  ssp_family = (keyword_set(miles))? 'MILES' : 'PEGASE.HR'
  sxaddpar,out_hdr_res,'SSPFAMIL',ssp_family

  sdssname = 'sdss'+basename
  out_file=outpath+dirsuff+'spSpec'+basename+'_results.fits'
  writefits,out_file,new_gal,h_new
  if(finite(chi2f[0]) eq 1) then begin
      res_tbl=[{ $
           SPEC_ID:0,       $
           NAME:sdssname,   $
           WAVE:wl_new,     $
           FLUX:galbin,     $
           ERROR:errbin,    $
           FIT:bestfit,     $
           FIT_COMP:allcompbestfitarr, $
           GOODPIXELS:flagGoodPix, $
           LSF_VEL:velbin,  $
           LSF_SIG:interpol(double(resData.DISPERSION)*velScaleSDSS,wlvac,wl_new), $
           LSF_SIG_REL:disbin,  $
           LSF_H3:disbin*0.0,   $
           LSF_H4:disbin*0.0,   $
           EM_LIN:emlines,  $
           V:velf[0,2*findgen(nlosvd)],     $
           E_V:velf[0,2*findgen(nlosvd)+1],   $
           SIG:disf[0,2*findgen(nlosvd)],   $
           E_SIG:disf[0,2*findgen(nlosvd)+1], $
           H3:h3f[0,2*findgen(nlosvd)],     $
           E_H3:h3f[0,2*findgen(nlosvd)+1],   $
           H4:h4f[0,2*findgen(nlosvd)],     $
           E_H4:h4f[0,2*findgen(nlosvd)+1],   $
           AGE:agef[0,0,*],   $
           E_AGE:agef[0,1,*], $
           MET:metf[0,0,*],   $
           E_MET:metf[0,1,*], $
           CHI2:chi2f[0],   $
           WEIGHT:ampf[0,0,*] $
      }]
;  res_tbl=[{ $
;           SPEC_ID:0,       $
;           NAME:sdssname,   $
;           WAVE:wl_new,     $
;           FLUX:galbin,     $
;           ERROR:errbin,    $
;           FIT:bestfit,     $
;           LSF_VEL:velbin,  $
;           LSF_SIG:disbin,  $
;           EM_LIN:emlines,  $
;           V:velf[0,0],     $
;           E_V:velf[0,1],   $
;           SIG:disf[0,0],   $
;           E_SIG:disf[0,1], $
;           H3:h3f[0,0],     $
;           E_H3:h3f[0,1],   $
;           H4:h4f[0,0],     $
;           E_H4:h4f[0,1],   $
;           AGE:agef[0,0],   $
;           E_AGE:agef[0,1], $
;           MET:metf[0,0],   $
;           E_MET:metf[0,1], $
;           CHI2:chi2f[0],   $
;           WEIGHT:ampf[0,0] $
;  }]
      mwrfits,res_tbl,out_file,out_hdr_res
  endif

  file_delete,tmpdir+'/'+tmpname,tmpdir+'/'+etmpname,tmpdir+'/'+'velmap.fts',$
              tmpdir+'/'+'dismap.fts',$
              tmpdir+'/'+'h3map.fts',tmpdir+'/'+'h4map.fts',tmpdir+'/'+'wl.fts',$
              tmpdir+'/'+'mags.fits',tmpdir+'/'+'errmags.fits',$
              /allow_nonexistent ;,/verb
endfor

if(deltmp eq 1) then spawn,'rmdir '+tmpdir

end