;VG

FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra

;+
; NAME:
;    dustem_mpfit_run
; PURPOSE:
;    function used to fit dustem model with dustem_mpfit_sed
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    res=dustem_mpfit_run(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
; INPUTS:
;    x         = necessary for interface but unused
;    params    = Dustem Model parameters normalized to initial values
; OPTIONAL INPUT PARAMETERS:
;    cc_sed    = SED color correction factors
;    _extra    = extra parameters for the plot routine
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    SED and model are plotted if !dustem_show_plot
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    res=dustem_mpfit_run(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
; MODIFICATION HISTORY:
;    Adapted to extinction and polarization by Vincent Guillet (IAS)
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

;note that x is not used.

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_mpfit_run'
  sed=0.
  goto,the_end
ENDIF
;stop

p_dim=p_min*(*(*!dustem_fit).param_init_values)

;CALCULATE THE CONTINUUM
IF !dustem_idl_continuum THEN cont=dustem_create_continuum()

;CALL THE DUSTEM_CREATE FUNCTIONS
;Actually, res is never used
;THIS line below was added by NF, but removed by JPB on Dec 5 2009
;dustem_set_params,p_dim         ;NF on October 2009 (params have to be updated before calling some functions !)
;dustem_create_functions,p_dim/(*!IDO),cont=cont,res=res
dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,res=res

;if ISRF has been modified and IONFRAC has not been modified, then call dustem_create_ionfrac
;NF on October 2009 : handling of DUSTEM_CREATE_IONFRAC has been changed
;this call only if IONFRACPAH master-keyword in GRAIN.DAT
;Removed for now (JPB). Use of ionfrac should be fully revisited.
IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()

;RUN THE MODEL AND GET THE SPECTRUM
st = dustem_run(p_dim)

;st_model=dustem_read_all(!dustem_soft_dir)

chi2_sed      = 0  
chi2_ext      = 0  
chi2_polext   = 0 
rchi2_sed     = 0
rchi2_ext     = 0
rchi2_polext  = 0
wrchi2_sed    = 0
wrchi2_ext    = 0
wrchi2_polext = 0
n_sed         = 0
n_ext         = 0
n_polext      = 0

tagnames=tag_names(!dustem_data)
dustem_out = [0]
for i=0,n_tags(!dustem_data)-1 do begin

    if !dustem_data.(i) eq ptr_new() then continue

    CASE tagnames(i) OF

	'SED'	: BEGIN

		
                dustem_sed = dustem_compute_sed(p_dim,st,cont=cont) 
                dustem_out = [ dustem_out, dustem_sed ]

		; Plot if needed
		IF !dustem_show_plot NE 0 THEN BEGIN

			window,1
  			;dustem_plot_fit_sed,st,dustem_sed,cont,_extra=extra,res=p_min*(*IDO),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
  			dustem_plot_fit_sed,st,dustem_sed,cont,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
			;  wait,0.5  ;Sorry for efficiency but needed to see the plot

	 	ENDIF

       		chi2_sed   = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
		n_sed = n_elements((*!dustem_data.sed).values)
       		rchi2_sed  = chi2_sed / n_sed
       		wrchi2_sed  = rchi2_sed * !fit_rchi2_weight.sed
       		print,"chi2 SED = ",chi2_sed,"  rchi2 SED = ",rchi2_sed

		END

	'EXT' 	: BEGIN 

		dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav) 
	 	dustem_out = [ dustem_out, dustem_ext ]

		; Plot if needed
		IF !dustem_show_plot NE 0 THEN BEGIN
		    	window,2
    			dustem_plot_ext,st,(*!dustem_params),yr=[1.e-6,1.e0],/ysty,xr=[1,400],/xsty
    			window,3
    			dustem_plot_ext,st,(*!dustem_params),/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty
	 	ENDIF

		dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
        	chi2_ext   = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2)
		n_ext = n_elements((*!dustem_data.ext).values)
		rchi2_ext  = chi2_ext / n_ext
        	wrchi2_ext  = rchi2_ext * !fit_rchi2_weight.ext
        	print,"chi2 EXT = ",chi2_ext,"  rchi2 EXT = ",rchi2_ext
		
	 	END

	'POLEXT': BEGIN

		dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav) 
		dustem_out = [ dustem_out, dustem_polext ]

		; Plot if needed
		IF !dustem_show_plot NE 0 THEN BEGIN

       			; Polarization fraction in the UV-IR
       			window,4
       			dustem_plot_polar,st,(*!dustem_params),yr=[0,0.08],/ysty,xr=[0.1,4],/xsty
			
       			; Polarisation curve (Serkowski)
       			window,5
       			dustem_plot_polar,st,(*!dustem_params),/UV,yr=[1e-3,3e-2],/ysty,xr=[0.1,10],/xsty

			; Also show prediction for polarized SED
			window,6
	        	tit='Predicted polarization fraction in emission'
        		ytit='Polarization fraction'
			xtit=textoidl('\lambda (\mum)')
			xrange = st.sed.wav
			Ngrains=n_tags(st.sed)-2
			colors=dustem_grains_colors(Ngrains)
	
        		plot_oi,xrange,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),/nodata,xtit=xtit,ytit=ytit,tit=tit,yr=[0,0.3],/xlog
        		for i = 0, Ngrains-1 do oplot,xrange,st.polsed.(i+1)/st.sed.(i+1),color=colors(i),psym=psym,linestyle=2
        		oplot,xrange,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),psym=psym,linestyle=2

	 	ENDIF

		dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
        	chi2_polext   = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
		n_polext = n_elements((*!dustem_data.polext).values)
        	rchi2_polext  = chi2_polext / n_polext
        	wrchi2_polext  = rchi2_polext * !fit_rchi2_weight.polext
        	print,"chi2 POLEXT = ",chi2_polext,"  rchi2 POLEXT = ",rchi2_polext

		END

  
        else :	stop

    ENDCASE

endfor

; Size distribution
;window,7
;fn = path+'out/SDIST.RES'
;OPENR, uu, fn, /get_lun
;tt = '#'
;WHILE STRPOS(tt,'#') EQ 0 do begin
;READF, uu, tt
;ENDWHILE
;ax = dblarr(ntype,nsz_max)
;ava = dblarr(ntype,nsz_max)
;FOR i=0,ntype-1 do begin
;READF,uu,tt
;for is=0,nsize(i)-1 do begin
;READF, uu, tx, ty
;ax(i,is) = tx
;ava(i,is) = ty
;endfor
;ENDFOR
;close,uu

;tit='Dust size distributions'
;yscl = 1.D29
;xr=[0.1,1e3]
;yr=[-2,3]
;xtit='a (nm)'
;ytit='Mass log( 10!u29!n n!dH!u-1!n a dm/da ) '
;plot_oi, ax(0,0:nsize(0)-1)*1.e7, alog10(ava(0,0:nsize(0)-1)*yscl), line=ls(i+1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit='DUSTEM'
;FOR i=1,ntype-1 DO begin
;	oplot, ax(i,0:nsize(i)-1)*1.e7, alog10(ava(i,0:nsize(i)-1)*yscl), line=ls(i+1)
;ENDFOR

print,"Total chi2 = ",chi2_sed+chi2_ext+chi2_polext
print,"Total weighted reduced chi2 driving mpfitfun = ", wrchi2_sed+wrchi2_ext+wrchi2_polext


; Suppress the artificial [0] at the beginning 
dustem_out = dustem_out[1:*]

;=== Caution, do not use message below here, even
;=== with the /info and /contine options, or the
;=== fit procedure will detect an error !
;=== use print instead.
IF defined(extra) THEN BEGIN
  IF keyword_set(extra.plot_it) THEN BEGIN
;    message,'params='+strtrim(p_min,2),/info,/continue
;    stop
    print,'dustem_mpfit_run: Current parameters values:',p_dim
  ENDIF
ENDIF

the_end:

RETURN,dustem_out

END