@pikaia.pro

; //////////////////////////////////////////////////////////////////////////////////////////
; Read the parameter file
; //////////////////////////////////////////////////////////////////////////////////////////

function read_param_file

  readcol, 'specfit_control.inp', names, values, format='(A,A)', delimiter='='

  param_name = [['Inner radius of the rim [AU]'],$
                ['Inner radius of the disk [AU]'],$
                ['Outer radius of the disk [AU]'],$
                ['Number of points in the rim'],$
                ['Number of points in the disk'],$
                ['Number of parameters to be fitted'],$
                ['Number of individuals in the pop.'],$
                ['Number of generations'],$
                ['Number of singificant digits'],$
                ['Crossover probability'],$
                ['Mutation mode'],$
                ['Initial mutation rate'],$
                ['Minimum mutation rate'],$
                ['Maximum mutation rate'],$
                ['Relative fitness differential'],$
                ['Reproduction plan'],$
                ['Elitism flag'],$
                ['Printed output (none/min/verbose)']]

   bound_name = [['Bounds for Trim'],$
                 ['Bounds for Tsurf_high'],$
                 ['Bounds for Tmid'],$
                 ['Bounds for qrim'],$
                 ['Bounds for qsurf'],$
                 ['Bounds for qmid'],$
                 ['Bounds for Tsurf_low']]


   param_val = fltarr(n_elements(param_name))*0.0 - 1
   bound_val = fltarr(2,n_elements(bound_name))*0.0 - 1

   for j = 0, n_elements(names)-1 do begin
      ; Scan the next line from the file (names[j] and values[j])
      idl = strpos(names[j], '<')
      idr = strpos(names[j], '>')
      text = strmid(names[j], idl+1, (idr-idl)-1)
      ; Check if it was a parameter
      for i=0, n_elements(param_name)-1 do begin
         if text eq param_name[i] then begin
            param_val[i] = float(values[j])
         endif
      endfor
      ; Check if it was a boundary
      for i=0, n_elements(bound_name)-1 do begin
         if text eq bound_name[i] then begin
            dum3 = strsplit(values[j], /extract)
            bound_val[0,i] = dum3[0]
            bound_val[1,i] = dum3[1]
         endif
      endfor
   endfor

   ctrl = param_val(6:17)
   rin  = param_val(0)
   rrim = param_val(1)
   rout = param_val(2)
   nr1  = param_val(3)
   nr2  = param_val(4)
   nvar = param_val(5)

   result = {ctrl:ctrl, nvar:nvar, rin:rin, rrim:rrim, rout:rout, nr1:nr1, nr2:nr2, pbound:bound_val}

   return, result

end
; //////////////////////////////////////////////////////////////////////////////////////////
; Initialize the fitting (read the opacities, PAH-templates, fitting parameter ranges, etc.)
; //////////////////////////////////////////////////////////////////////////////////////////
pro read_all, source_list_file=source_list_file, nr=nr

common param_bounds, pbound
common opacities, kappa, ndust, qval_fname, pah_fname, npah, pah_prof
common obs_params, dist, rstar, tstar, nwav, owav, oflux, ouflux, fit_range, source_name, rsub
common grid_params, r1_min, r1_max, r2_min, r2_max, nr1, nr2
common directories, model_dir

model_dir        = './BestFit/'
spectra_dir      = '' ;'./SPECTRA/'
thisname         = 'read_all'

@natconst.pro
;
; Read the spectral name list
;
readcol, source_list_file, dum1, dum2, dum3, dum4, dum5, dum6, dum7, dum8, dum9, dum10, $
         format='(I,A,A,A,F,F,F,F,F,F)'

ii = where(dum1 eq nr)
if ii[0] ge 0 then begin
   if n_elements(ii) gt 1 then begin
      print, 'ERROR! Multiple files have been found with the same Nr'
      stop
   endif else begin
      spec_fname   = dum2(ii(0))
      source_name  = dum3(ii(0))
      sed_fname    = dum4(ii(0))
      fit_range    = [dum5(ii(0)), dum6(ii(0))]
      rstar        = dum7(ii(0)) * rs
      tstar        = dum8(ii(0))
      dist         = dum9(ii(0)) * pc
      rsub         = dum10(ii(0))
  endelse
endif else begin
   print, 'ERROR! No spectrum with the specified Nr has been found! '
   stop
endelse

;
; Read the spectrum itself, select the interval [fit_range[0], fit_range[1]]
;
readcol, spectra_dir+spec_fname, wav, flux, uflux, skipline=1
ii = where(wav ge fit_range[0] and wav le fit_range[1])
nwav   = n_elements(ii)
if nwav le 1 then begin
  print, thisname, ': Not enough datapoints in fitrange'
  stop
endif
owav   = wav(ii)
oflux  = flux(ii)
ouflux = uflux(ii)
;
; Read the opacity and PAH template file names
;
  openr, 1, 'opacity_files.inp'
  qval_path = ''
  readf, 1, qval_path
  readf, 1, ndust, format='(I)'
  qval_fname = ''
  for i=0, ndust-1 do begin
     dum = ''
     readf, 1, dum
     ; Each entry consists of a number (1/0) and the filename
     ; The file will only be used if the number != 0
     sdum = strsplit(dum, /extract)
     dum1 = fix(sdum[0])
     dum2 = sdum[1]
     if dum1 eq 1 then qval_fname = [qval_fname, dum2]
  endfor
  qval_fname = qval_fname(1:n_elements(qval_fname)-1)
  ndust = n_elements(qval_fname)
;
; Skip any number of empty lines
;
  dum1 = ''
  dum2 = ''
  while dum1 eq '' do readf, 1, dum1
;
; Read the PAH template filenames, dum1 is their directory
;
  readf, 1, dum2
  pah_path = dum1
  npah = fix(dum2)
  pah_fname = ''
  for i=0, npah-1 do begin
     dum = ''
     readf, 1, dum
     ; Each entry consists of a number (1/0) and a filename
     ; The fiel is only used if number != 0
     sdum = strsplit(dum, /extract)
     dum1 = fix(sdum[0])
     dum2 = sdum[1]
     if dum1 eq 1 then pah_fname = [pah_fname, dum2]
  endfor
  if n_elements(pah_fname) eq 1 and pah_fname(0) eq '' then begin
     npah = 0
     pah_fname = ''
  endif else begin
     pah_fname = pah_fname(1:n_elements(pah_fname)-1)
     npah = n_elements(pah_fname)
  endelse
  close, 1
  free_lun, 1

;
; Compatibility check with the Fortran90 version
;
  ii = strpos(qval_path, "'")
  if ii(0) ge 0 then begin
     qval_path = strsplit(qval_path, "'", /extract)
     if n_elements(qval_path) ne 1 then qval_path = qval_path(1)
  endif

  ii = strpos(pah_path, "'")
  if ii(0) ge 0 then begin
     pah_path = strsplit(pah_path, "'", /extract)
     if n_elements(pah_path) ne 1 then pah_path = pah_path(1)
  endif

;
; Read the opacities
;
  kappa = fltarr(nwav, ndust)
  for iq=0, n_elements(qval_fname)-1 do begin

     readcol, qval_path[0] + qval_fname[iq], dum1, dum2, skipline=1, /silent
     readcol, qval_path[0] + qval_fname[iq], nw, gsize, gdens, numline=1, /silent
     replace = where(dum2 lt 0)
     if n_elements(replace) gt 0 and replace[0] ne -1 then dum2[replace] = 0.0
     ; Interpolate the Qext to the wavelength grid of the data
     idum2 = interpol(dum2, dum1, owav)
     ; Compute mass absorption coefficients from the Qext
     kappa(*,iq) = 3.0 * idum2 / (4. * gsize[0] * 1d-4 * gdens[0])
     print, qval_fname[iq], ' .. reading .... DONE', format='(A40, A30)'
  endfor

;
; Read the PAH-templates
;
  if npah gt 0 then begin
     pah_prof = fltarr(nwav, npah)
     for ip=0, n_elements(pah_fname)-1 do begin
        readcol, pah_path[0] + pah_fname[ip], dum1, dum2, skipline=1, /silent, format='(F,F)'
        replace = where(dum2 lt 0)
        if n_elements(replace) gt 0 and replace[0] ne -1 then dum2[replace] = 0.0
        ; Interpolate the PAH profiles to the wavelength grid of the data
        pah_prof[*,ip] = interpol(dum2, dum1, owav)
        print, pah_fname[ip], ' .. reading .... DONE', format='(A40, A30)'
     endfor
  endif

end

; //////////////////////////////////////////////////////////////////////////////////////////
; Write the results into a file
; //////////////////////////////////////////////////////////////////////////////////////////
pro write_output, ci, comps, spar, redchi, logfile_name

common opacities, kappa, ndust, qval_fname, pah_fname, npah, pah_prof
common obs_params, dist, rstar, tstar, nwav, owav, oflux, ouflux, fit_range, source_name
common grid_params, r1_min, r1_max, r2_min, r2_max, nr1, nr2

@natconst

;   openw, 1, 'fitresults.log', /append
  openw, 1, logfile_name, /append

  norm_high = total(ci(3:ndust+2))
  norm_low = total(ci(ndust+3:2*ndust+2))


    printf, 1
    printf, 1,  ' *********************************************** '
    printf, 1,  '    Input parameters :                           '
    printf, 1,  ' *********************************************** '
    printf, 1,  '       Name of the star            : ', source_name
    printf, 1,  '       Radius of the star          : ', Rstar/RS
    printf, 1,  '       Temperature of the star     : ', Tstar
    printf, 1,  '       Distance                    : ', dist/pc
    printf, 1,  '       Fitted wavelength range     : ', fit_range(0), ' ',fit_range(1)

    printf, 1,  '       Inner disk radius           : ', r1_min/AU
    printf, 1,  '       Radius of the rim           : ', r1_max/AU
    printf, 1,  '       Outer disk radius           : ', r2_max/AU

    printf, 1,  ' *********************************************** '
    printf, 1,  '    The fitted parameters : '
    printf, 1,  ' *********************************************** '
    printf, 1,  '       Temperature of the rim      : ', spar(0)
    printf, 1,  '       Temperature of the surface  : ', spar(1)
    ;printf, 1,  '       Temperature of the surface  : ', spar(6)
    printf, 1,  '       Temperature or the midplane : ', spar(2)

    printf, 1,  '       Prim                        : ', spar(3)
    printf, 1,  '       Psurf                       : ', spar(4)
    printf, 1,  '       Pmidpl                      : ', spar(5)
    printf, 1,  ' *********************************************** '
    printf, 1,  '    Reduced chisquare  : ', redchi
    printf, 1,  ' *********************************************** '
    printf, 1,  '       Contribution of the Star     : ', ci(0)
    printf, 1,  '       Contribution of the Rim      : ', ci(1)
    printf, 1,  '       Contribution of the Midplane : ', ci(2)
    printf, 1,  '       Contribution of the Surface (high)  : ', norm_high
    printf, 1,  '       Contribution of the Surface (low)   : ', norm_low
    printf, 1,  ' *********************************************** '
    printf, 1,  '    Surface (high):'
    printf, 1,  ' *********************************************** '

    for i=0, ndust-1 do begin
       if norm_high eq 0.0 then mf = 0.0 $
       else mf = (ci(i+3)/norm_high*100.) ;round(ci(i+3)/norm_high*100.)
       printf, 1, qval_fname(i), mf, format='(A40,"  ", F8.4)'
    endfor
    printf, 1,  ' *********************************************** '
    printf, 1,  '    Surface (low):'
    printf, 1,  ' *********************************************** '
    for i=0, ndust-1 do begin
       if norm_low eq 0. then mf = 0.0 $
       else mf = (ci(i+3+ndust)/norm_low*100.) ;round(ci(i+3+ndust)/norm_low*100.)
       printf, 1, qval_fname(i), mf, format='(A40,"  ", F8.4)'
    endfor
    printf, 1,  ' *********************************************** '
    printf, 1,  '     PAH:'
    printf, 1,  ' *********************************************** '
    if npah gt 0 then begin
       for ip=0, npah-1 do begin
          printf, 1, pah_fname(ip), ci(ip+3+2*ndust), format='(A40,"  ", F)'
       endfor
    endif

  close, 1
end

;==============================================================================
; Write an ascii table of the measured spectrum and the various components of
; the best fit model. The columns in the table are as follows: wavelength [mu],
; observed flux [Jy], total model flux [Jy], continuum in model [Jy], PAH [Jy],
; dust [Jy]. Two lines of header, commented with '#' will precede the table. In
; the first line will be the source name, while the second is for the column
; headers.
; Input:  model, total model flux
;         dust, dust emission from surface layers
;         pah, PAH emission
;         cont, summed continuum (i.e. midplane, star, ...)
; Return: 0 on success, 1 on error
;------------------------------------------------------------------------------

; function writeResults, model, dust, pah, cont
function writeResults, model, dust, cont

COMMON obs_params, dist, rstar, tstar, nwav, owav, oflux, ouflux, fit_range, $
                   source_name
COMMON directories, model_dir

fname  = model_dir + '/' + source_name + '_results.dat'
; chead  = 'wav [mu] flux [Jy] model [Jy]  cont. [Jy]  PAH [Jy] dust [Jy]'
; format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4)"
chead  = 'wav [mu] flux [Jy] model [Jy]  cont. [Jy]  dust [Jy]'
format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4)"


openw, fd, fname, /get_lun, width=200
printf, fd, '#  ' + source_name
printf, fd, '# ' + chead
for i = 0, nwav - 1 do begin
  ; printf,fd, owav[i], oflux[i], model[i], cont[i], pah[i], dust[i], format=format
  printf,fd, owav[i], oflux[i], model[i], cont[i], dust[i], format=format
endfor
close, fd
free_lun, fd

return, 0
end

; //////////////////////////////////////////////////////////////////////////////////////////
; Planck function
; //////////////////////////////////////////////////////////////////////////////////////////
function planck, wav, temp

  cc     = 2.9979d10
  kk     = 1.3807d-16
  hh     = 6.6262d-27
  ss     = 5.6703d-5
  au     = 1.496d13

  nu     = 2.9979d14/wav

  planck = (2.d*hh*nu^3.d/cc^2) / (exp(hh*nu/(kk*temp))-1)

  return, planck
end

; //////////////////////////////////////////////////////////////////////////////////////////
;   Calculate the source function (integrating over radius)
; //////////////////////////////////////////////////////////////////////////////////////////

function get_source_func_rad, wav=wav, spar=spar

common obs_params, dist, rstar, tstar, nwav, owav, oflux, ouflux, fit_range, source_name
common grid_params, r1_min, r1_max, r2_min, r2_max, nr1, nr2

@natconst.pro

trim0  = spar(0)
tsurf0_high = spar(1)
tmid0  = spar(2)
prim   = spar(3)
psurf  = spar(4)
pmid   = spar(5)
;tsurf0_low = spar(6) ;JV 2019!!!

conv   = 1d0 / dist^2. * 1d23

;
; Create the spatial grid and the corresponding temperatures
;

r1 = r1_min * (r1_max/r1_min)^(dindgen(nr1)/(nr1-1.))
r2 = r2_min * (r2_max/r2_min)^(dindgen(nr2)/(nr2-1.))


;
; Get the boundary between the low and high temp surf component
;

;rbound = r2_min * (tsurf0_low / tsurf0_high)^(1./psurf(0))
rbound = r2_max ;JV 2019!!!
;
; Compute the temperature grid(s)
;

iih = where(r2 le rbound)
iil = where(r2 gt rbound)
nr2l = n_elements(iil)
nr2h = n_elements(iih)

trim  = trim0[0]  * (r1/r1_min)^prim[0]
tsurf_high = tsurf0_high[0] * (r2/r2_min)^psurf[0]
;tsurf_low = tsurf0_low[0] * (r2/r2_min)^psurf[0] ;JV 2019!!!
tmid  = tmid0[0]  * (r2/r2_min)^pmid[0]

;
; Calculate the emission from the star
;

starsurf = !dpi * rstar^2.
sflux = planck(wav, tstar) * starsurf * conv

;
; Calculate the source function for each annulus
;

brim       = dblarr(nr1,nwav)
bsurf_high = dblarr(nr2h,nwav)
;bsurf_low  = dblarr(nr2l,nwav)
bmid       = dblarr(nr2,nwav)

; Inner rim
for ir=0, nr1-2 do begin
   da   = !dpi * (r1(ir+1)*r1(ir+1) - r1(ir)*r1(ir))
   temp = (trim(ir) + trim(ir+1))*0.5
   brim(ir+1, *) = planck(wav, temp(0)) * da * conv
endfor

; High temperature atmosphere
if(nr2h gt 1) then begin
    for ir=0, nr2h-2 do begin
        da = !dpi * (r2[ir+1]^2. - r2[ir]^2.)
        temp = (tsurf_high[ir] + tsurf_high[ir+1])*0.5
        bsurf_high[ir+1, *] = planck(wav, temp[0]) * da * conv
    endfor
endif

; Low temperature atmosphere  ;JV 2019!!!
;if (nr2 gt nr2h) then begin
;    for ir=nr2h-1, nr2-2 do begin
;        da = !dpi * (r2(ir+1)^2. - r2(ir)^2.)
;        temp = (tsurf_low(ir) + tsurf_low(ir+1))*0.5
;        bsurf_low(ir+1-nr2h, *) = planck(wav, temp(0)) * da * conv
;    endfor
;endif

; Midplane
for ir=0, nr2-2 do begin
   da = !dpi * (r2(ir+1)^2. - r2(ir)^2.)
   temp = (tmid(ir) + tmid(ir+1))*0.5
   bmid(ir+1, *) = planck(wav, temp(0)) * da * conv
endfor

;
; Sum up the contribution of each annulus to each wavelength
;

for inu=0, nwav-1 do begin
   ; The rim
   brim(0,inu)  = total(brim(1:nr1-1,inu), 1)
   buffer = size(bsurf_high)
   ; The high temperature atmosphere
   if (nr2h gt 1) then bsurf_high(0,inu) = total(bsurf_high(1:nr2h-1,inu), 1) $
    else bsurf_high[0, inu] = 0.d
   ; The low temperature atmosphere
   ;if (nr2 - nr2h gt 1) then bsurf_low(0,inu) = total(bsurf_low(1:nr2-nr2h-1,inu), 1) $ ;JV 2019!!!
    ;else bsurf_low[0, inu] = 0.d
   ; The midplane
   bmid(0,inu)  = total(bmid(1:nr2-1,inu), 1)
endfor

;
; Create a structure to return
;

res = {wav:wav, bstar:sflux, brim:brim, bsurf_high:bsurf_high, bsurf_low:bsurf_high*0.0, bmid:bmid}  ;JV 2019!!!

return, res
end


; //////////////////////////////////////////////////////////////////////////////////////////
;  Scale the paramters from [0,1] to [par_min, par_max]
; //////////////////////////////////////////////////////////////////////////////////////////
function scale_params, par

common param_bounds, pbound


spar = fltarr(n_elements(par))
for i=0, n_elements(par)-1 do begin
   spar(i) = pbound(0,i) + (pbound(1,i) - pbound(0,i)) * par(i)
;   spar[i] = (par[i] - pbound[0,i])/(pbound[1,i] - pbound[0,i])
endfor

return, spar

end


; //////////////////////////////////////////////////////////////////////////////////////////
;  Calculate the model
; //////////////////////////////////////////////////////////////////////////////////////////

function get_model, par, ci=ci, comps=comps

common opacities, kappa, ndust, qval_fname, pah_fname, npah, pah_prof
common obs_params, dist, rstar, tstar, nwav, owav, oflux, ouflux, fit_range, source_name

;
; Scale the parameters
;

spar = scale_params(par)

;
; Calculate the corresponding source function
;

src = get_source_func_rad(wav=owav, spar=spar)

;
; Calculate the emissivities
;

fsurf_high = dblarr(nwav, ndust)
fsurf_low = dblarr(nwav, ndust)
for idust=0, ndust-1 do begin
   fsurf_high(*,idust) = kappa(*,idust) * src.bsurf_high(0,*)
   ;fsurf_low(*,idust) = kappa(*,idust) * src.bsurf_low(0,*)
endfor

;
; Prepare for the linear decomposition
;

bq = dblarr(nwav, 2*ndust+npah+3)
bq(*,0) = src.bstar
bq(*,1) = src.brim(0,*)
bq(*,2) = src.bmid(0,*)
bq(*,3:ndust+2) = fsurf_high
bq(*,ndust+3:2*ndust+2) = fsurf_low
if npah gt 0 then bq(*,2*ndust+3:2*ndust+npah+2) = pah_prof

bound = dblarr(2, 2*ndust+npah+3)
bound(0,*) = 0d0
bound(1,*) = 1d4

bound(0,0) = 0.8
bound(1,0) = 1.2

bound(0,1) = 0.
bound(1,1) = 1.

bound(0,2) = 0.
bound(1,2) = 1.



;
; Do the linear decomposition
;

bk_bq = bq
bk_flux = oflux

bvls, bk_bq, bk_flux, bound, ci, iter=iter, ierr=ierr

;
; Calculate the model function
;

fmod = bq # ci

;
;  Create a structure to return
;

comps = bq
res = {fmod:fmod, ci:ci}

return, res

end

; //////////////////////////////////////////////////////////////////////////////////////////
;  Calculate the chisquare
; //////////////////////////////////////////////////////////////////////////////////////////

function get_chi, n, par

common obs_params, dist, rstar, tstar, nwav, owav, oflux, ouflux, source_name
common opacities, kappa, ndust, qval_fname, pah_fname, npah

res = get_model(par)

chi = total((oflux - res.fmod)^2. / ouflux^2.)
; We fit 2*ndust + npah + 3 parameters ci and the BB temperatures of rim, midplane
; and the two atmosphere components
dof = nwav - (2*ndust + npah + 3 + 4)

rchi = chi / dof

return, 1./rchi

end

; //////////////////////////////////////////////////////////////////////////////////////////
;  PIKAIA driver routine
; //////////////////////////////////////////////////////////////////////////////////////////

function specfit_tdga, nr=nr, source_list_file=source_list_file, logfile_name=logfile_name, rsub=rsub

common param_bounds, pbound
common opacities, kappa, ndust, qval_fname, pah_fname, npah, pah_prof
common obs_params, dist, rstar, tstar, nwav, owav, oflux, ouflux, fit_range, source_name
common grid_params, r1_min, r1_max, r2_min, r2_max, nr1, nr2
COMMON directories, model_dir

@natconst.pro

if not keyword_set(nr) then nr=1
if not keyword_set(source_list_file) then source_list_file='source_list_fit_fulap.inp'
if not keyword_set(logfile_name) then logfile_name = 'fitresults.log'
;
; Initialize the fitting reading everything
;
read_all, nr=nr, source_list_file=source_list_file

;
; Read the parameter input file
;

par = read_param_file()

pbound = par.pbound

; rin    = par.rin
; rrim   = par.rrim
rin    = rsub
rrim   = rsub*1.1
rout   = par.rout
nr1    = par.nr1
nr2    = par.nr2
nvar   = par.nvar
ctrl   = par.ctrl

;
; Spectrum we want to fit
;

r1_min         = rin * au
r1_max         = rrim * au
r2_min         = rrim * au
r2_max         = rout * au

;--------------------------------------------------------------------------------------------
; RUN THE GENETIC FITTING ROUTINE TO OBTAIN THE BEST TEMPERATURE DISTRIBUTIONS
;--------------------------------------------------------------------------------------------
pikaia, nvar, 'get_chi', ctrl, bestpar, fitness, status
model = get_model(bestpar, ci=ci, comps=comps)
spar  = scale_params(bestpar)

chi = total((oflux - model.fmod)^2. / ouflux^2.)
dof = nwav - (2*ndust + 3 + 4)
rchi = chi / dof
;--------------------------------------------------------------------------------------------
; WRITE THE RESULTS INTO A FILE
;--------------------------------------------------------------------------------------------
write_output, ci, comps, spar, rchi, logfile_name
;--------------------------------------------------------------------------------------------
; PLOT THE RESULTS
;--------------------------------------------------------------------------------------------

fcomps = fltarr(nwav,2*ndust+3+npah)
for i=0, 2*ndust+2+npah do fcomps(*,i) = comps(*,i) * ci(i)
sflux_high = total(fcomps(*,3:ndust+2), 2)
;sflux_low  = total(fcomps(*,ndust+3:2*ndust+2), 2)
cont       = total(fcomps(*,0:2), 2)
; pahs       = total(fcomps(*,2*ndust+3:2*ndust+2+npah),2)

norm_high = total(ci(3:ndust+2))
norm_low = total(ci(ndust+3:2*ndust+2))
mfrac_high = ci(3:ndust+2)/norm_high*100.
mfrac_low  = ci(ndust+3:2*ndust+2)/norm_low*100.

set_plot, 'ps'
!P.FONT = 0

r = [0,1,1,0,0, 1.0,1.0,0.5,0.0,0.5,0.0, 0.0,0.5,0.5, 0.75,0.25,0.25,  0.75,0.50,0.50, 0.75,0.25,0.25] * 255
g = [0,1,0,1,0, 0.5,0.0,1.0,1.0,0.0,0.5, 0.5,0.0,0.5, 0.25,0.75,0.25,  0.50,0.75,0.25, 0.25,0.75,0.50] * 255
b = [0,1,0,0,1, 0.0,0.5,0.0,0.5,1.0,1.0, 0.5,0.5,0.0, 0.25,0.25,0.75,  0.25,0.25,0.75, 0.50,0.50,0.75] * 255
tvlct, r,g,b

plotname = './results_plot/' + source_name + '.eps'
device, filename=plotname, ysize=18.0, xsize=24.0,yoffset=24, $
        xoffset=0, /color, /landscape, /encapsulated

plot, owav, oflux, psym=10, title=source_name
oplot, owav, model.fmod, color=2
oplot, owav, cont,       color=3
oplot, owav, sflux_high, color=4
;oplot, owav, sflux_low,  color=6
for j=0,ndust-1 do BEGIN
  oplot, owav, fcomps(*,3+j+2), color=6+j
endfor

; oplot, owav, pahs,       color=7
; al_legend, ['model', 'cont.', 's_high', 's_low', 'pah'], linestyle=[0,0,0,0,0], color=[2,3,4,6,7]
al_legend, ['model', 'cont.', 'surf. dust'], linestyle=[0,0,0], color=[2,3,4]

device, /close
set_plot, 'X'

cmd = 'eps2png -f -r 90 '+plotname
print,cmd
spawn,cmd,cmd_output
print,cmd_output

; dump = writeResults(model.fmod, sflux_high + sflux_low, pahs, cont)
;dump = writeResults(model.fmod, sflux_high + sflux_low,  cont)
;write results
fname  = model_dir + '/' + source_name + '_results.dat'
; chead  = 'wav [mu] flux [Jy] model [Jy]  cont. [Jy]  PAH [Jy] dust [Jy]'
; format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4)"
chead  = ' wav_mu    flux_Jy      model_Jy     cont_Jy      dust_Jy       '
format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ',$)"

openw, fd, fname, /get_lun, width=200
printf, fd, '#  ' + source_name
printf, fd, '# ' + chead, format="(A,$)"
for j=0,ndust-1 do BEGIN
  printf,fd,qval_fname[j],format="(A,' ',$)"
endfor
printf,fd,''
for i = 0, nwav - 1 do begin
  ; printf,fd, owav[i], oflux[i], model[i], cont[i], pah[i], dust[i], format=format
  printf,fd, owav[i], oflux[i], model.fmod[i], cont[i], sflux_high[i], format=format
  for j=0,ndust-1 do BEGIN
    printf,fd, fcomps(i,3+j+2), format="(E12.4, ' ',$)"
  endfor
  printf,fd,''
endfor
close, fd
free_lun, fd

end