;+
; NAME:            ULY_TGM_EVAL
; 
; PURPOSE:           
;                  Evaluate a TGM model spectrum.
; 
; USAGE:
;                  array = uly_tgm_eval(eval_data, para)
; 
; ARGUMENTS:
;   EVAL_DATA:     The interpolated stellar spectrum model coefficients, 
;                  it contains three dimensions of data, specifying to
;                  warm, hot and cold stellar regions respectively, in
;                  fact, they are the coefficients of the evaluating 
;                  polynomial function.
;
;   PARA:          Input 3 atmospheric parameters (log(Teff), Logg, [Fe/H]),
;                  they will be passed to ULY_TGM_MODEL_PARAM to
;                  compute out the parameters which will be applied
;                  to the polynomial function for evaluating
;                  the model spectrum at this given stellar
;                  atmospheric parameters (Teff, log g and [Fe/H]).
;
; DESCRIPTION:
;      Evaluate a TGM model spectrum which is interpolated as a
;      function of atmospheric parameters (effective temperature,
;      surface gravity and [Fe/H]). Each wavelength point is evaluated
;      by using polynomials.
;
;      This function is automatically called by the minimization routines:
;      ULY_FIT, ULY_FIT_LIN to evaluate a TGM model (see ULY_TGM).
; 
; OUTPUT:          
;                  Evaluated stellar model spectrum               
;
; REQUIRED FUNCTION: 
;                  ULY_TGM_MODEL_PARAM, compute the parameters 
;       
; AUTHOR:
;                  Yue WU, 2008/06/18
;
;-
; CATEGORY:        ULY_TGM
;------------------------------------------------------------------------------
function uly_tgm_model_param, version, teff, gravi, fehi  
                                                                 
compile_opt idl2, hidden
on_error, 2

if version le 2 then param = dblarr(23,3) $
else if version eq 3 then param = dblarr(26,3) $
else message, 'Version'+string(version)+' of TGM not supported'

for i = 0,2 do begin
    teffc = teff*1d0                            
    if (version eq 2 or version eq 3) and i eq 2 then teffc = teff + 0.1
    grav = gravi * 1d0
    feh = fehi * 1d0

    tt = teffc/0.2 
    tt2 = tt^2-1. 
                                                                 
    param[0,i] = 1.                                                      
    param[1,i] = tt                     
    param[2,i] = feh                                                     
    param[3,i] = grav                                                    
    param[4,i] = tt^2                   
    param[5,i] = tt*(tt2)               
    param[6,i] = tt2*tt2                
    param[7,i] = tt*feh                                                  
    param[8,i] = tt*grav                                                 
    param[9,i] = tt2*grav                                                
    param[10,i] = tt2*feh                                                
    param[11,i] = grav^2                                                 
    param[12,i] = feh^2                                                  
    param[13,i] = tt*tt2^2              
    param[14,i] = tt*grav^2                                              
    param[15,i] = grav^3                                                 
    param[16,i] = feh^3                                                  
    param[17,i] = tt*feh^2                                               
    param[18,i] = grav*feh                                               
    param[19,i] = grav^2*feh                                             
    param[20,i] = grav*feh^2                                             
    if version eq 1 then begin
       param[21,i] = exp(tt) 
       param[22,i] = exp(tt*2) 
    endif else if version eq 2 then begin
       param[21,i] = exp(tt) - 1d - tt*(1d + tt/2d + tt^2/6d + tt^3/24d + tt^4/120d)
       param[22,i] = exp(tt*2) - 1d - 2D*tt*(1d + tt + 2d/3d*tt^2 + tt^3/3d + tt^4*2d/15d)

    endif else if version eq 3 then begin
       param[21,i] = tt*tt2*grav            
       param[22,i] = tt2*tt2*grav            
       param[23,i] = tt2*tt*feh            
       param[24,i] = tt2*grav^2
       param[25,i] = tt2*grav^3
    endif else message, 'unsupported version of TGM interpolator '+string(version)

endfor

return, param                                                    
end   

;------------------------------------------------------------------------------
function uly_tgm_eval, eval_data, para

spec_coef = eval_data.spec_coef
spec_npix = (size(spec_coef, /DIM))[0]

teff = alog10(exp(para[0]))-3.7617                        
grav = para[1]-4.44   

param = uly_tgm_model_param(eval_data.version, teff, grav, para[2])
np =  (size(param,/DIM))[0] - 1

if (teff le alog10(9000.)-3.7617) then t1 = spec_coef[*,0:np,0] # param[*,0]
if (teff ge alog10(7000.)-3.7617) then t2 = spec_coef[*,0:np,1] # param[*,1]
if (teff le alog10(4550.)-3.7617) then t3 = spec_coef[*,0:np,2] # param[*,2]


if (teff le alog10(7000.)-3.7617) then begin
    if (teff gt alog10(4550.)-3.7617) then begin 
        tgm_model_evalhc = t1
    endif else if (teff gt alog10(4000.)-3.7617) then begin 
        q = (teff-alog10(4000.)+3.7617)/(alog10(4550.)-alog10(4000.))
        tgm_model_evalhc = q*t1+(1.-q)*t3
    endif else begin
        tgm_model_evalhc = t3
    endelse
endif else if (teff ge alog10(9000.)-3.7617) then begin
    tgm_model_evalhc = t2
endif else begin
    q = (teff-alog10(7000.)+3.7617)/(alog10(9000.)-alog10(7000.))
    tgm_model_evalhc = q*t2+(1.-q)*t1
endelse


; Resample in wavelength if necessary (when option REBIN_COEF is not given)
if eval_data.sampling ne eval_data.mod_samp or $
   eval_data.start ne eval_data.mod_start or $
   eval_data.step ne eval_data.mod_step or $
   eval_data.npix ne spec_npix then begin
   spec = uly_spect_alloc(DATA=tgm_model_evalhc, START=eval_data.mod_start, $
                          STEP=eval_data.mod_step, SAMPLING=eval_data.mod_samp)
   wrange = eval_data.start + [0d, eval_data.npix * eval_data.step]
   if eval_data.sampling eq 1 then wrange = exp(wrange)
   c = 299792.458d              ; Speed of light in km/s
   velscale = eval_data.step * c
   if eval_data.sampling eq 1 then $
      spec = uly_spect_logrebin(spec, velscale, WAVERANGE=wrange, /OVER) $
   else $
      spec = uly_spect_linrebin(spec, eval_data.step, WAVERANGE=wrange, /OVER)
   tgm_model_evalhc = *spec.data
   uly_spect_free, spec
endif

if eval_data.calibration eq 'C' then begin ; multiply by a black-body spectrum
  n = n_elements(tgm_model_evalhc)
  if eval_data.sampling eq 1 then $
    wavelength = exp(eval_data.start + dindgen(n) * eval_data.step) $
  else $
    wavelength = eval_data.start + dindgen(n) * eval_data.step 
  w = 5550.
  c3 = 1.43883d8 / 5550d0 / exp(para[0])
  c1 = 3.74185d19 / 5550d0^5
  if c3 lt 50. then bbm = (c1 / (exp(c3)-1.)) else bbm = (c1 * exp(-c3))

  c3 = 1.43883d8 / wavelength / exp(para[0])
  c1 = 3.74185d19 / wavelength^5 / bbm

  n1 = where(c3 lt 50, cnt, COMPLEM=n2)
  if cnt gt 0 then tgm_model_evalhc[n1] *= (c1[n1] / (exp(c3[n1])-1.))
  if cnt lt n then tgm_model_evalhc[n2] *= (c1[n2] * exp(-c3[n2]))

endif

; Convolve LOSVD in case giving lsf_file
; (if REBIN_COEF is given, the LSF injection shall be made at initialization.)
if n_elements(eval_data.lsf) gt 0 and eval_data.lsf ne 'no_lsf' then begin
    spec = uly_spect_alloc(DATA=tgm_model_evalhc, START=eval_data.start, $
                           STEP=eval_data.step, SAMPLING=1)
    uly_spect_lsfconvol, eval_data.lsf, spec
    tgm_model_evalhc = *spec.data
    uly_spect_free, spec
endif

return, tgm_model_evalhc

end