FUNCTION compute_TLS_poly_gaussian_fit, temp, wave, parameter_set=parameter_set

; To get informations on this function please read Paradis et al. (2014). 
; wave: wavelength(s) (in microns)
; temp: temperature(s) (in K)
;parameter_set='diffuse' or 'compact_source' or 'standard'
; CAUTION : This function returns Qabs(temp,wave) in arbitrary units!!
; dustem_planck_function.pro is available on the Dustem wrapper webpage:
; http://dustemwrap.irap.omp.eu/ by downloading the DustEM wrapper

;Written by Deborah Paradis. 
; If you have any questions, please contact me: dparadis@irap.omp.eu


   coefs=make_array(5,5)
   aa=make_array(6)
IF parameter_set EQ 'diffuse' then begin
   coefs(*,0)=[3.33830, -5.36544e-3, 1.65381e-6, 1.63754e-11, -6.88816e-12]
   coefs(*,1)=[-2.98801e-4, 1.49526e-5, -1.42524e-7, 1.41845e-12, 6.72034e-13]
   coefs(*,2)=[8.54787e-6, -2.17671e-7, 1.49001e-9, 1.03925e-11, -2.48183e-14]
   coefs(*,3)=[-3.25327e-8, 1.06390e-9, 4.20971e-12, -1.37459e-13, 2.39831e-16]
   coefs(*,4)=[4.55497e-11, -1.10996e-12, -3.84120e-14, 4.36257e-16, -6.77351e-19]
   aa(*,0)=[-0.00050, 0.07585, 5.36111, 100.19785, 13.47090, 499.96309]
ENDIF

IF parameter_set EQ 'compact_source' then begin
   coefs(*,0)=[3.34100, -6.57687e-3, 1.67055e-5, -7.42107e-8, 7.97846e-11]
   coefs(*,1)=[-1.78243e-3, 8.93472e-5, -1.09494e-6, 3.82566e-9, -2.98224e-12]
   coefs(*,2)=[4.27571e-5, -1.54310e-6, 1.94520e-8, -5.15802e-11, 2.51914e-14]
   coefs(*,3)=[-2.53887e-7, 1.06665e-8, -1.19930e-10, 2.44969e-13, -2.12459e-17]
   coefs(*,4)=[5.47140e-10, -2.43342e-11, 2.47310e-13, -3.67727e-16, -2.21749e-19]
   aa(*,0)=[-0.00126, 0.14795, 4.52068, 91.28719, 9.50025, 451.84591]
ENDIF

IF parameter_set EQ 'standard' then begin
   coefs(*,0)=[3.33042, -5.49209e-3, 2.14115e-6, -2.14173e-9, -7.90430e-12]
   coefs(*,1)=[-2.84206e-4, 1.12572e-5, -1.11099e-7, -7.48509e-11, 9.32649e-13]
   coefs(*,2)=[9.23806e-6, -1.28278e-7, 1.17408e-9, 1.19308e-11, -3.01610e-14]
   coefs(*,3)=[-3.20750e-8, 4.97171e-10, 5.92135e-12, -1.52179e-13, 2.83618e-16]
   coefs(*,4)=[3.66206e-11, 1.88461e-13, -4.33218e-14, 4.82933e-16, -7.97139e-19]
   aa(*,0)=[-0.00062, 0.09271, 5.14453, 90.71464, 12.36827, 484.87299]
ENDIF

Nt=n_elements(temp)
Nw=n_elements(wave)
logqabs=dblarr(Nt,Nw)
qabs=logqabs
V=qabs
ggauss=V
tmin=4.30
dt=6.46369e-1
wmin=2.00076
dw=2.70167e-3
w=alog10(wave)
t=temp

for j=0,Nt-1 do begin
   for i=0,Nw-1 do begin
      logqabs(j,i)=1.d*(coefs(0,0)+coefs(0,1)*((t(j)-tmin)/dt)+coefs(0,2)*((t(j)-tmin)/dt)^2+coefs(0,3)*((t(j)-tmin)/dt)^3+coefs(0,4)*((t(j)-tmin)/dt)^4  $
                        +coefs(1,0)*((w(i)-wmin)/dw)+coefs(1,1)*((w(i)-wmin)/dw)*((t(j)-tmin)/dt)+coefs(1,2)*((w(i)-wmin)/dw)*((t(j)-tmin)/dt)^2+coefs(1,3)*((w(i)-wmin)/dw)*((t(j)-tmin)/dt)^3+coefs(1,4)*((w(i)-wmin)/dw)*((t(j)-tmin)/dt)^4 $
                        +coefs(2,0)*((w(i)-wmin)/dw)^2+coefs(2,1)*((w(i)-wmin)/dw)^2*((t(j)-tmin)/dt)+coefs(2,2)*((w(i)-wmin)/dw)^2*((t(j)-tmin)/dt)^2+coefs(2,3)*((w(i)-wmin)/dw)^2*((t(j)-tmin)/dt)^3+coefs(2,4)*((w(i)-wmin)/dw)^2*((t(j)-tmin)/dt)^4 $
                        +coefs(3,0)*((w(i)-wmin)/dw)^3+coefs(3,1)*((w(i)-wmin)/dw)^3*((t(j)-tmin)/dt)+coefs(3,2)*((w(i)-wmin)/dw)^3*((t(j)-tmin)/dt)^2+coefs(3,3)*((w(i)-wmin)/dw)^3*((t(j)-tmin)/dt)^3+coefs(3,4)*((w(i)-wmin)/dw)^3*((t(j)-tmin)/dt)^4 $
                        +coefs(4,0)*((w(i)-wmin)/dw)^4+coefs(4,1)*((w(i)-wmin)/dw)^4*((t(j)-tmin)/dt)+coefs(4,2)*((w(i)-wmin)/dw)^4*((t(j)-tmin)/dt)^2+coefs(4,3)*((w(i)-wmin)/dw)^4*((t(j)-tmin)/dt)^3+coefs(4,4)*((w(i)-wmin)/dw)^4*((t(j)-tmin)/dt)^4)
      V(j,i)=((t(j)-aa(4))/aa(2))^2 + ((w(i)-aa(5))/aa(3))^2
      ggauss(j,i)=aa(0) + aa(1) * EXP(-V(j,i)/2.) 
;stop
      qabs(j,i)=10.^((logqabs(j,i))-ggauss(j,i))*dustem_planck_function(t(j),10.^(w(i)))
   endfor
endfor

return, qabs
END