FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=_extra

;+
; NAME:

;    dustem_mpfit_run_vg

; PURPOSE:
;    function used to fit dustem model with dustem_mpfit_sed
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:

;    res=dustem_mpfit_run_vg(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_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])

;  ODIFICATION 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

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

IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()

;RUN THE MODEL AND GET THE SPECTRUM

dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values)

st = dustem_run(p_dim)

st_model=dustem_read_all(!dustem_soft_dir)

chi2_sed      = 0  
chi2_ext      = 0  
chi2_polext   = 0 

;chi2_polsed   = 0 
;chi2_polfrac   = 0

chi2_qsed      =0
chi2_used      =0 
chi2_qext      =0
chi2_uext      =0 

rchi2_sed     = 0
rchi2_ext     = 0
rchi2_polext  = 0

;rchi2_polsed  = 0
;rchi2_polfrac  = 0

rchi2_qsed     =0
rchi2_used     =0
rchi2_qext     =0
rchi2_uext     =0

wrchi2_sed    = 0
wrchi2_ext    = 0
wrchi2_polext = 0

;wrchi2_polsed = 0
;wrchi2_polfrac = 0

wrchi2_qsed    =0
wrich2_used    =0
wrchi2_qext    =0
wrich2_uext    =0

n_sed         = 0
n_ext         = 0
n_polext      = 0

;n_polsed      = 0
;n_polfrac      = 0

n_qsed         =0
n_used         =0
n_qext         =0
n_uext         =0

win = 0

tagnames=tag_names(!dustem_data)
dustem_out = [0]

ind_data_tag=dustem_data_tag(count=count_data_tag)

;if exist('./noplot') then !dustem_show_plot = 0 else if !dustem_show_plot eq 0 then !dustem_show_plot = 1

;stop

FOR i_tag=0,count_data_tag-1 DO BEGIN

;stop

;print,i_tag,ind_data_tag(i_tag)
;print,tagnames(ind_data_tag(i_tag))

    CASE tagnames(ind_data_tag(i_tag)) OF

		'SED'	: BEGIN

	        dustem_sed = dustem_compute_sed(p_dim,st)

	        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
	       	message,"chi2 SED = "+strtrim(chi2_sed,2)+"  rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed
			;print,'wl ',(*!dustem_data.sed).wav
			;print,"Delta SED per wl: ",(*!dustem_data.sed).values-dustem_sed
			;print,"error per wl: ",(*!dustem_data.sed).sigma
			;print,"chi2 per wl: ",((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2
	        alpha = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $
	              		/ total(((*!dustem_data.sed).values)^2/(*!dustem_data.sed).sigma^2)
	        beta = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $
	             		/ total(dustem_sed^2/(*!dustem_data.sed).sigma^2)

	        print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta

	        (*!dustem_fit).chi2=chi2_sed
	        (*!dustem_fit).rchi2=rchi2_sed

			; Plot if needed
	        IF !dustem_show_plot NE 0 THEN BEGIN

	            ;window,win,title='DUSTEM WRAP (SED)' 

	            ;cgwindow,WMULTI=[2,2,1],/CURRENT

                
                
                dustemwrap_plot,p_dim,st,dustem_sed,dustem_polsed,dustem_qsed,dustem_polfrac,_extra=_extra
                ;dustemwrap_plott,p_dim,st,dustem_sed,dustem_polsed,dustem_qsed,dustem_polfrac,_extra=_extra
	            ;dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2,fpol=fpol

		 	ENDIF

	        dustem_out = [ dustem_out, dustem_sed ]

		END

		'EXT' 	: BEGIN 

    		
    		;ADDING PLUGIN TO SPECTRUM----------------

    		IF ptr_valid(*!dustem_plugin) THEN BEGIN

    		for i=0L,n_tags(*!dustem_scope)-1 do begin

            if total(strsplit(*(*!dustem_scope).(i),'+',/extract) eq 'ADD_EXT') then st.ext.ext_tot+=(*(*!dustem_plugin).(i))

            endfor

            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;," weighted rchi2 EXT=",wrchi2_ext
	        	;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2

			; Plot if needed
			IF !dustem_show_plot NE 0 THEN BEGIN

				win+=1
	                	;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win
				xr=[0.3,40]
				yr=[1.e-26,1.e-21]
				dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
				dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2

				win+=1
	                	;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
				xr=[0.2,10]
				yr=[0,3e-21]
				dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
				dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2

			ENDIF
			

			;dustem_plot_ext,st,/print_Rv ;commented this recently

		 	dustem_out = [ dustem_out, dustem_ext ]

		 	END

		'POLEXT': BEGIN

			
			;---------------the procedure within this block has not been updated for a vert long time

			IF !dustem_show_plot NE 0 and !run_circ then begin
	                   win+=1
	                 
				dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03]
			endif

			;--------------------------------------------------------------------------------------
			
			

			IF !run_lin then begin

                    ;compute_polext needs updating. It's been months. 

                    dustem_polext = dustem_compute_polext(p_dim,st,dustem_qext,dustem_uext,Q_ext,U_ext)

	        		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;," weighted rchi2 POLEXT=",wrchi2_polext
	        		;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2

	        		; Fit in (pmax,lmax,K)
	        		x=st.polext.wav
	        		logx=alog(x)
	        		logy=alog(st.polext.ext_tot)

	        		; Whittet+1992 (Table 1)
	        		Uband = 0.36
	        		Bband = 0.43
	        		Vband = 0.55
	        		Rband = 0.63
	        		Iband = 0.78
	        		Jband = 1.21
	        		Hband = 1.64
	        		Kband = 2.04
	        		bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband]
	        		nbands = n_elements(bands)
	        		jband = indgen(nbands)*0
	        		for k=0,nbands-1 do begin
	                		xx=min(abs(st.polext.wav-bands(k)),j)
	                		jband(k) = j
	        		endfor
	        		;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
	        		res = poly_fit(logx(jband),logy(jband),2)
			
	        		print
	        		print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)')
	        		print,"-------------------------"
	        		K    = -res(2)
	        		lmax = exp(res(1)/(2*K))
	        		pmax = exp(res(0)+K*alog(lmax)^2)
	        		print,"* lmax = ",lmax," micron"
	        		print,"* K    = ",K
	        		print,"* pmax = ",pmax," for NH=1e21"
	        		print
			
				; Plot if needed
				IF !dustem_show_plot NE 0 THEN BEGIN

;  	       			; Polarization fraction in the UV-IR
;  	                win+=1
; 	                           
;  	       			dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,50],/xsty,win=win;multi=0,ps=filename
; ; 					

	   				; Polarisation curve (Serkowski)

 					win+=1
 					;xr=[0.3,1e4]
 					xr=[1e-4,3.5]
 					yr=[1e-30,5e-24]

 					
 					dustem_plot_polext,st,p_dim,dustem_polext,aligned=st_model.align.grains.aligned,win=win
 					
;  					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,/donotclose,multi=1
;  					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,multi=2

;  					win+=1
;  					dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
;  					

	                                
	            ENDIF

				;dustem_plot_polar,st,/print_ratio,cont=cont
	                     
				dustem_out = [ dustem_out, dustem_polext ]

			ENDIF

			END

; 		'POLFRAC': BEGIN ;Deprecate this. Fitting POLFRAC does not make any sense. It is a biased quantity. 
;     		if !run_lin then begin
;                 dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here.      
;                 ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat        
;                 ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths    
;                 
;                 ;dustem_sed and dustem_polsed might not have the same length thus this block.
;                 
;                 If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin                            
;                     if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin 
;                         dustem_polsed_x = dustem_polsed
;                         dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
;                         ;matching the filter data points
;                         IF count_filt ne 0 then begin
;                             for i=0L,count_filt-1 do begin    
;                                 j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
;                                 if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0))                                                 
;                             endfor                                        
;                         endif
;                         ;matching the spectrum data points    
;                         IF count_fspec ne 0 then begin
;                             for i=0L,count_fspec-1 do begin
;                                 j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
;                                 if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0))                                    
;                             endfor
;                         endif     
;                     endif ELSE begin
;                         
;                         dustem_polsed_x = dustem_sed
;                         dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
;                         ;matching the filter data points
;                         IF count_filt ne 0 then begin
;                             for i=0L,count_filt-1 do begin    
;                                 j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
;                                 if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0))                                                 
;                             endfor                                        
;                         endif
;                         ;matching the spectrum data points    
;                         IF count_fspec ne 0 then begin
;                             for i=0L,count_fspec-1 do begin
;                                 j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
;                                 if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0))                                    
;                             endfor
;                         endif
;                     
;                     ENDELSE
;                 endif ELSE BEGIN
;                     dustem_polsed_x = dustem_polsed 
;                     dustem_sed_x = dustem_sed
;                 ENDELSE
;                 
;                 dustem_polfrac = dustem_polsed_x/dustem_sed_x 
; 	        	chi2_polfrac   = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
; 				n_polfrac = n_elements((*!dustem_data.polfrac).values)
; 	        	rchi2_polfrac  = chi2_polfrac / n_polfrac
; 	        	wrchi2_polfrac  = rchi2_polfrac * !fit_rchi2_weight.polfrac
; 	        	message,"chi2 POLFRAC = "+strtrim(chi2_polfrac,2)+"  rchi2 POLFRAC = "+strtrim(rchi2_polfrac,2),/continue    ;," weighted rchi2 POLFRAC=",wrchi2_polfrac
;                 dustem_out = [ dustem_out, dustem_polfrac ]
;                 ;fpol=dustem_polfrac ; problem because extra structure does not get updated so you actually need the 'extra' procedure you haven't finihed coding.
;             endif
;OLD DEPRECATED  VERSION.             

;             
; 			IF !run_lin then begin
; 
; 	            dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron)
; 				dustem_sed_tmp = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
; 				dustem_sed = dustem_polsed
; 				ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
; 				for i = 0, count_filt-1 do begin
; 					j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
; 					if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
; 				endfor 
; 				dustem_polfrac = dustem_polsed/dustem_sed
; 
; 			; We normalize at 353GHz
; 			;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
; 			;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)
; 
; 	        	chi2_polfrac   = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
; 				n_polfrac = n_elements((*!dustem_data.polfrac).values)
; 	        	rchi2_polfrac  = chi2_polfrac / n_polfrac
; 	        	wrchi2_polfrac  = rchi2_polfrac * !fit_rchi2_weight.polfrac
; 	        	print,"chi2 POLFRAC = ",chi2_polfrac,"  rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
; 	        	;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2
; 
; 	            dustem_out = [ dustem_out, dustem_polfrac ]
; 
; 			; Plot if needed
; 				IF !dustem_show_plot NE 0 THEN BEGIN
; 
; 				; Also show prediction for polarized P/I
; 					win+=1
; 	       			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,yr=[0,0.25],/ysty,xr=[10,5d3];dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,yr=[0,0.25],/ysty,xr=[10,5d3]
; 
; 		 		ENDIF
; 
; 		 	ENDIF
; 
; 			END


;-----------------------------------------------------------------

; 		'POLSED': BEGIN
; 		    
; 		    ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
; 			
; 			IF !run_lin THEN BEGIN
;                 ;make polsed only used for when polsed is fitted. Otherwise it's the fitting of Q_SED AND USED that fits the polarization data
; 	        	dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
; 	        
; 	        	chi2_polsed   = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
; 				n_polsed = n_elements((*!dustem_data.polsed).values)
; 	        	rchi2_polsed  = chi2_polsed / n_polsed
; 	        	wrchi2_polsed  = rchi2_polsed * !fit_rchi2_weight.polsed
; 	        	message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+"  rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
; 	        	print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
; 	        	;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
;               
; 	            dustem_out = [ dustem_out, dustem_polsed ] ; classic command
; 		        ;-----------------------------------------------
; 	         
; 				; Plot if needed
; 				IF !dustem_show_plot NE 0 THEN BEGIN
; 					win+=1
; ; 					
; 					xr=[10,8e3]
; 					yr=[1e-34,1e-28]

; 					

; 					dustem_plot_polsed,st,p_dim,dustem_polsed,aligned=st_model.align.grains.aligned,win=win
; ; 					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1
; ; 					dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
; 
;                     					
; 					
; 				ENDIF
; 			ENDIF
; 		END

        'QSED': BEGIN 

		        toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) 

		        ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
		        dustem_qsed(-1)=dustem_qsed(-2)
			    stop

			    chi2_qsed   =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2) 
			    n_qsed = n_elements((*!dustem_data.qsed).values)
			    rchi2_qsed  = chi2_qsed / n_qsed
	        	wrchi2_qsed  = rchi2_qsed * !fit_rchi2_weight.qsed 
	            dustem_out = [ dustem_out, dustem_qsed ] ; classic command
	            message,"chi2 QSED = "+strtrim(chi2_qsed,2)+"  rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
	        	print,"chi2 per wl: ",((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2
		        
		      
		        
		        if !dustem_show_plot NE 0 then begin

		        fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
    	        fact1 = 1.00E-20;(3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting   

; 	            fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17

;               fact1=1.E-20

; 	            

		        win+=1
                window  ,win ,title='DUSTEM WRAP (QSED)'
		        titstq='Stokes Q Polarization Intensity'

                ytitstq=textoidl('Q_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')

                xtit=textoidl('\lambda (\mum)')

                xr=[10,1e6]
               

                ;normpol1=(*!dustem_data.qsed).values*fact1

                yr=[1e-28,1e-17]
                ;yr=[1e-34,1e-23]

                indq=where(Q_spec lt 0, countq)

                ;indq=where((*!dustem_data.qsed).values lt 0, countq)
                

                ylog=1
                if countq ne 0 then begin

;                 ylog=0
;                 yr=[-1E-23,-1e-19];e-31 
                    Q_spec=-Q_spec
                    dustem_qsed =-dustem_qsed
                    ((*!dustem_data.qsed).values)=-((*!dustem_data.qsed).values)

                endif

                normpol=Q_spec*fact
                normpol1=dustem_qsed*fact1

                cgplot,st.polsed.wav,Q_spec*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)'
                cgoplot,st.polsed.wav,Q_spec*fact,color='black'

                cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1,color=cgColor('Tomato')
                cgoplot,(*!dustem_data.qsed).wav,dustem_qsed*fact1,psym=6,color='red',symsize=2
                

                cgplot,st.polsed.wav,Q_spec*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,st.polsed.wav,Q_spec*fact/normpol,color='black'

                cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato')
                        
                endif

		    END   

        'USED': BEGIN
		    

                toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)

                stop
 			    chi2_used   =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2) 

 			    n_used = n_elements((*!dustem_data.used).values)
 			    rchi2_used  = chi2_used / n_used

  	        	wrchi2_used  = rchi2_used * !fit_rchi2_weight.used

 	            dustem_out = [ dustem_out, dustem_used] ; classic command
 	            message,"chi2 USED = "+strtrim(chi2_used,2)+"  rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed

  	        	print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2
  
 	            
 	            if !dustem_show_plot ne 0 then begin
 	            fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
     	        fact1 =1.00E-20;(3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting   

 	            win+=1

                window  ,win ,title='DUSTEM WRAP (USED)'                 
                titstu='Stokes U Polarization Intensity'

                ytitstu=textoidl('U_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')

                xtit=textoidl('\lambda (\mum)')

     	        xr=[10,1e6]
     	        normpol2=U_spec*fact

                normpol3=dustem_used*fact1
                ;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10]	

                  ;conditions here are on computed spectra and they should be on the data itself.
                  ;This even impairs the visualization of the data. Correct this!!!!!
                  
                  yr=[1e-28,1e-17];[1e-34,1e-23]
                  indu=where(U_spec lt 0, countu)
                  ;indu=where((*!dustem_data.used).values lt 0, countu)
                  ylog=1
                  if countu ne 0 then begin
                  ylog=0
                  yr=[-2e-31,5e-31]  
                  endif 

                cgplot,st.polsed.wav,U_spec*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
                cgoplot,st.polsed.wav,U_spec*fact,color='black'

                cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,color='Teal',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1,color=cgColor('Teal')
                cgoplot,(*!dustem_data.used).wav,dustem_used*fact1,psym=6,color='red',symsize=2
                
                
                

                cgplot,st.polsed.wav,U_spec*fact/normpol2,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,st.polsed.wav,U_spec*fact/normpol2,color='black'

                cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,color='Teal',psym=16,symsize=1,thick=2
                err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1/normpol3,color=cgColor('Teal')
                                
                ;cgwindow, 'plot', 
                	 

 	            endif
  		END
;  		
;  		

 		
 		
  		'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
                   
 			    chi2_qext   =total(((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2) 
   			    n_qext = n_elements((*!dustem_data.qext).values)
   			    rchi2_qext  = chi2_qext / n_qext
  	        	wrchi2_qext  = rchi2_qext * !fit_rchi2_weight.qext
   	            dustem_out = [ dustem_out, dustem_qext] ; classic command
   	            message,"chi2 QEXT = "+strtrim(chi2_qext,2)+"  rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  	        	print,"chi2 per wl: ",((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2
                if !dustem_show_plot ne 0 then begin
	            
	            win+=1
                window  ,win ,title='DUSTEM WRAP (QEXT)'	                 
                titqex='Q-polarized Cross-section'
                ytitqex=textoidl('\sigma_{Q}/N_H (cm^2/H)')
                xtit=textoidl('\lambda^{-1} (\mum^{-1})')
    	        ;xr=[0.3,1e4]
    	        xr=[1e-4,3.5]
    	          
                yr=[1e-34,5e-24]
                indqex=where(Q_ext lt 0, countqex)
                ylog=1
                if countqex ne 0 then begin
                ylog=0
                yr=[-9e-27,5e-27]  
                endif 
                normpol = st.polext.ext_tot * 0. + 1.d21  
                cgplot,1/st.polext.wav,Q_ext,xtit='',ytit=ytitqex,tit=titqex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
                cgoplot,1/st.polext.wav,Q_ext/normpol,color='black'
                cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,color='Indian Red',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,yrms=3.*((*!dustem_data.qext).sigma)/2/normpol,color=cgColor('Indian Red')
                cgoplot,1/(*!dustem_data.qext).wav,dustem_qext/normpol,psym=6,color='red',symsize=2
                
                
                
                cgplot,1/st.polsed.wav,Q_ext/Q_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,1/st.polsed.wav,Q_ext/Q_ext,color='black'
                cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,color='Indian Red',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,yrms=3.*((*!dustem_data.qext).sigma)/2.,color=cgColor('Indian Red')
                    
            	            endif                  
     	           
  		END
;  		
  		'UEXT': BEGIN
  
 			    
 			    chi2_uext   =total(((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2) 
   			    n_uext = n_elements((*!dustem_data.uext).values)
   			    rchi2_uext  = chi2_uext / n_uext
  	        	wrchi2_uext  = rchi2_uext * !fit_rchi2_weight.uext
   	            dustem_out = [ dustem_out, dustem_uext] ; classic command
   	            message,"chi2 UEXT = "+strtrim(chi2_uext,2)+"  rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  	        	print,"chi2 per wl: ",((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2
                if !dustem_show_plot ne 0 then begin
	            
	            win+=1
                window  ,win ,title='DUSTEM WRAP (UEXT)'	                 
                tituex='U-polarized Cross-section'
                ytituex=textoidl('\sigma_{U}/N_H (cm^2/H)')
                xtit=textoidl('\lambda^{-1} (\mum^{-1})')
    	        ;xr=[0.3,1e4]
    	        xr=[1e-4,3.5]
    	          
                yr=[1e-34,5e-24]
                induex=where(U_ext lt 0, countuex)
                ylog=1
                if countuex ne 0 then begin
                ylog=0
                yr=[-9e-27,5e-27]  
                endif 
                normpol = st.polext.ext_tot * 0. + 1.d21  
                cgplot,1/st.polext.wav,U_ext,xtit='',ytit=ytituex,tit=tituex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
                cgoplot,1/st.polext.wav,U_ext/normpol,color='black'
                cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,color='Steel Blue',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
                cgoplot,1/(*!dustem_data.uext).wav,dustem_uext/normpol,psym=6,color='red',symsize=2
                
                
                
                cgplot,1/st.polsed.wav,U_ext/U_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1    
                cgoplot,1/st.polsed.wav,U_ext/U_ext,color='black'
                cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,color='Steel Blue',psym=16,symsize=1,thick=2
                err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
                    
            	            endif                  
     	        
                
                  
 	           
  		END
 		

	    ELSE : IF tagnames(ind_data_tag(i_tag)) NE 'POLFRAC' THEN stop 

    ENDCASE

ENDFOR

DOF = n_elements(p_min)
total_chi2 = 0
total_rchi2 = 0
total_n = 0
if !fit_rchi2_weight.ext ne 0 then begin
	total_n += n_ext
	total_chi2 +=  chi2_ext
	total_rchi2 += rchi2_ext
endif
if !fit_rchi2_weight.sed ne 0 then begin
	total_n += n_sed
	total_chi2 += chi2_sed
	total_rchi2 += rchi2_sed
endif
if !run_pol then begin
if !fit_rchi2_weight.polext ne 0 then begin
	total_n += n_polext
	total_chi2 += chi2_polext
	total_rchi2 += rchi2_polext
endif

; if !fit_rchi2_weight.polsed ne 0 then begin
; 	total_n += n_polsed 
; 	total_chi2 += chi2_polsed 
; 	total_rchi2 += rchi2_polsed 
; endif
; if !fit_rchi2_weight.polfrac ne 0 then begin
; 	total_n += n_polfrac
; 	total_chi2 += chi2_polfrac
; 	total_rchi2 += rchi2_polfrac
; endif

if !fit_rchi2_weight.qsed ne 0 then begin
	total_n += n_qsed
	total_chi2 += chi2_qsed
	total_rchi2 += rchi2_qsed
endif
if !fit_rchi2_weight.used ne 0 then begin
	total_n += n_used
	total_chi2 += chi2_used
	total_rchi2 += rchi2_used

endif

if !fit_rchi2_weight.qext ne 0 then begin
	total_n += n_qext
	total_chi2 += chi2_qext
	total_rchi2 += rchi2_qext
endif
if !fit_rchi2_weight.uext ne 0 then begin
	total_n += n_uext
	total_chi2 += chi2_uext
	total_rchi2 += rchi2_uext
endif

endif

total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $
            + chi2_sed * !fit_rchi2_weight.sed 
if !run_pol then begin
	total_wchi2 	+= chi2_polext  * !fit_rchi2_weight.polext $

;              		 + chi2_polsed  * !fit_rchi2_weight.polsed $
;              		 + chi2_polfrac * !fit_rchi2_weight.polfrac $

             		 + chi2_qsed * !fit_rchi2_weight.qsed $
             		 + chi2_used * !fit_rchi2_weight.used $
             		 + chi2_qext * !fit_rchi2_weight.qext $
             		 + chi2_uext * !fit_rchi2_weight.uext  

endif

message,'======================================================================',/continue

print,"Total weighted chi2 driving mpfitfun = ", total_wchi2
print,"Nb of data points = ", total_n
print,"DOF = ", DOF

print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed,rchi2_polext;,rchi2_polsed,rchi2_polfrac
print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed,wrchi2_polext;,wrchi2_polsed,wrchi2_polfrac
if !run_pol then 	print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!fit_rchi2_weight.polext $;,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $ 

else 			print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed
print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)

print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original

message,'======================================================================',/continue

;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF)
; Reduced chi2 is like if there were the same nb of data points for each curve
;print,"Mean reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac) * total_n/(total_n-DOF)

; Plot size distribution
win+=1
xr=[0.3,5e3]
yr=[0.1,100]

;IF !dustem_show_plot NE 0 THEN dustem_plot_sdist,xr=xr,yr=yr,ps=filename,win=win;,/donotclose

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

;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_vg: Current parameters values:',p_dim
;  ENDIF
;ENDIF

the_end:

RETURN,dustem_out

END