update from deborah

Jean-Michel Glorian
1 parent 427f1205
Showing 22 changed files with 939 additions and 276 deletions Show diff stats
src/idl/dustem_cc.pro
src/idl/dustem_compute_polsed.pro
src/idl/dustem_data_tag.pro
src/idl/dustem_grains_colors.pro
src/idl/dustem_init.pro
src/idl/dustem_mpfit_data.pro
src/idl/dustem_mpfit_run.pro
src/idl/dustem_plot_ext.pro
src/idl/dustem_plot_fpol.pro
src/idl/dustem_plot_nuinu_em.pro
src/idl/dustem_plot_polar.pro
src/idl/dustem_plot_sdist.pro
src/idl/dustem_read_align.pro
src/idl/dustem_read_all_res.pro
src/idl/dustem_set_data.pro
src/idl/mpfit.pro
src/idl/mpfitfun.pro
src/idl_misc/init_healpix.pro
src/idl_misc/isarray.pro
src/idl_misc/plotsym.pro
@@ -78,8 +78,8 @@ instrus=dustem_filter2instru(filter_names)
 inst_tags=tag_names((*!dustem_filters))
  
 ;DEBUG VG
-j=where(instrus eq 'HFI',count)
-if count gt 0 then instrus(j) = 'HFI_V3P02'
+;j=where(instrus eq 'HFI',count)
+;if count gt 0 then instrus(j) = 'HFI_V3P02'
  
 Nbands=n_elements(filter_names)
 sed=dblarr(Nbands)
@@ -53,6 +53,7 @@ spec = stp.em_tot * fact
 if not isarray(spec) THEN stop 
 ;stop
 ;COMPUTE THE MODEL SED
+
 dustem_polsed = (*!dustem_data.polsed).values * 0.
  
 ;stop
@@ -4,7 +4,8 @@ ntags=n_tags(!dustem_data)
 mm=lonarr(ntags)
  
 FOR i=0L,ntags-1 DO BEGIN
-  IF ptr_valid(!dustem_data.(i)) THEN mm(i)=1
+   IF ptr_valid(!dustem_data.(i)) THEN mm(i)=1
+
 ENDFOR
 ind=where(mm EQ 1,count)
  
-FUNCTION dustem_grains_colors,Ngrains
+FUNCTION dustem_grains_colors,Ngrains,ls,ps=ps
  
-colmax=210 & colmin=30
-colors=findgen(Ngrains)/(1.*Ngrains-1)*(colmax-colmin)+colmin
+fixed_colors = [0,210,50,210,50,30,50,80,120,150,180,210,230,250]
+;last_color = n_elements(fixed_colors)-1
+;colors = [0,fixed_colors(last_color-Ngrains+1:last_color)]
+colors = [0,fixed_colors(1:Ngrains)]
+
+ls = [ [ 0,4,3,5,1 ], [ 0,4,3,5,1 ],[ 0,4,3,5,1 ] ]
+
+if keyword_set(ps) then begin
+	; Postscript
+       	!y.thick = 5
+       	!x.thick = 5
+       	!p.thick = 7
+       	!p.charsize = 1.3
+       	!p.charthick = 5
+       	!x.ticklen = 0.04
+endif else begin
+	j=where(colors eq 0)
+	colors(j)=255
+endelse
  
 RETURN,colors
  
@@ -81,7 +81,8 @@ if keyword_set(pol) then begin
                                 polsed: ptr_new(),      $
                                 polfrac: ptr_new()       $
                          }
-	rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.}
+        rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.}
+    
 endif else begin
 	defsysv, '!dustem_data', {    $ ;Data to fit
                                 sed:    ptr_new(),      $
@@ -293,7 +294,7 @@ endif ; /noerase
 IF keyword_set(wraptest) THEN BEGIN
   dir_out=!dustem_dat
   dustem_write_all,st_model,dir_out
-  stop
+ ; stop
   st=dustem_run()
   IF keyword_set(plot_it) THEN BEGIN
     loadct,13
@@ -88,13 +88,14 @@ for j=0,n_elements(fields)-1 do begin
 		endelse
 		;print,"WEIGHT",i,j,weight,n_elements((*!dustem_data.(i)).sigma),"/",n_elements(wav)
  
-		instruc += ', (*!dustem_data.(' + strtrim(i,2) + ')).'+fields(j)+'*'+strtrim(weight,2)
+		instruc += ', (*!dustem_data.(' + strtrim(ind_data_tag(i),2) + ')).'+fields(j)+'*'+strtrim(weight,2)
 		;print,"INSTRUC",instruc
  
 	endfor
 	instruc += '])(1:*)'	; from 1 to exclude first 0
 	;print,'-> INSTRUC: ',instruc
-	toto=execute(instruc)	; Variables wav, values & sigma are created
+        toto=execute(instruc)	; Variables wav, values & sigma are created
+
      endfor
  
 p_min = dustem_mpfitfun(func_name,wav,values,sigma, $
-;VG
-
 FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra
  
 ;+
 ; NAME:
-;    dustem_mpfit_run
+;    dustem_mpfit_run_vg
 ; 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])
+;    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
@@ -32,7 +30,7 @@ FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra
 ; PROCEDURE:
 ;    None
 ; EXAMPLES
-;    res=dustem_mpfit_run(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help])
+;    res=dustem_mpfit_run_vg(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
@@ -42,16 +40,17 @@ FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra
  
 ;note that x is not used.
  
+;COMMON MPFIT_ERROR, ERROR_CODE
+
 IF keyword_set(help) THEN BEGIN
-  doc_library,'dustem_mpfit_run'
+  doc_library,'dustem_mpfit_run_vg'
   sed=0.
   goto,the_end
 ENDIF
-;stop
  
 p_dim=p_min*(*(*!dustem_fit).param_init_values)
+;Run dustem with as an interface to mpfitfun
  
-;CALCULATE THE CONTINUUM
 IF !dustem_idl_continuum THEN cont=dustem_create_continuum()
  
 ;CALL THE DUSTEM_CREATE FUNCTIONS
@@ -66,173 +65,361 @@ dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,res=
 ;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)
+st_model=dustem_read_all(!dustem_soft_dir)
  
 chi2_sed      = 0  
 chi2_ext      = 0  
 chi2_polext   = 0 
+chi2_polsed   = 0 
+chi2_polfrac   = 0 
 rchi2_sed     = 0
 rchi2_ext     = 0
 rchi2_polext  = 0
+rchi2_polsed  = 0
+rchi2_polfrac  = 0
 wrchi2_sed    = 0
 wrchi2_ext    = 0
 wrchi2_polext = 0
+wrchi2_polsed = 0
+wrchi2_polfrac = 0
 n_sed         = 0
 n_ext         = 0
 n_polext      = 0
+n_polsed      = 0
+n_polfrac      = 0
+
+win = 0
  
 tagnames=tag_names(!dustem_data)
 dustem_out = [0]
-for i=0,n_tags(!dustem_data)-1 do begin
+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
  
-    if !dustem_data.(i) eq ptr_new() then continue
+for i_tag=0,count_data_tag-1 do begin
  
-    CASE tagnames(i) OF
+;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,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
+		; We normalize to 353
+		;x=min(((*!dustem_data.sed).wav-850)^2,j353)
+		;dustem_sed = dustem_sed/dustem_sed(j353)*(*!dustem_data.sed).values(j353)
  
        		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
+       		print,"chi2 SED = ",chi2_sed,"  rchi2 SED = ",rchi2_sed;," 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_plot_nuinu_em,st,dustem_sed,cont,/print_radiance,p_dim=p_dim
+               
+		; Plot if needed
+                IF !dustem_show_plot NE 0 THEN BEGIN
+                   window,0
+                   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
+    			;ytit='Total intensity'
+                        ;xtit=textoidl('\lambda (\mum)')
+			;win+=1
+			;xr=[10,5e3]
+			;yr=[5e-33,5.e-28]
+			;nodiff = 0
+			;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1,nodiff=nodiff
+			;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2,nodiff=nodiff
+
+	 	ENDIF
+
+                dustem_out = [ dustem_out, dustem_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
+        	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
+
+	 	dustem_out = [ dustem_out, dustem_ext ]
+
 	 	END
  
 	'POLEXT': BEGIN
  
-		dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav) 
-		dustem_out = [ dustem_out, dustem_polext ]
+		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
+
+			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;," 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)
+        		;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
+
+        		; 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,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,4],/xsty,win=win
+				
+   				; Polarisation curve (Serkowski)
+				win+=1
+				xr=[0.1,50]
+				yr=[5e-25,5e-23]
+				dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,/donotclose,multi=1
+				dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,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
+
+		IF !run_lin then begin
+
+                dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont) 
+		dustem_sed_tmp = dustem_compute_sed(p_dim,st,cont=cont) 
+		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
  
-       			; 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
+			; Also show prediction for polarized P/I
+			win+=1
+       			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,yr=[0,0.25],/ysty,xr=[10,5d3]
  
 	 	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
+	 	ENDIF
+
+		END
+
+	'POLSED': BEGIN
+
+		IF !run_lin then begin
+
+                dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont) 
+
+		; We normalize to 353
+		;x=min(((*!dustem_data.polsed).wav-850)^2,j353)
+		;dustem_polsed = dustem_polsed/dustem_polsed(j353)*(*!dustem_data.polsed).values(j353)
+
+        	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
+        	print,"chi2 POLSED = ",chi2_polsed,"  rchi2 POLSED = ",rchi2_polsed;," 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 ]
+
+		; Plot if needed
+		IF !dustem_show_plot NE 0 THEN BEGIN
+
+			win+=1
+			xr=[10,5e3]
+			yr=[1e-34,1e-28]
+			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,xr=xr,yr=yr,/xsty,/donotclose,multi=1
+			dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,xr=xr,yr=yr,/xsty,multi=2
+
+		ENDIF
+
+		ENDIF
  
 		END
  
-  
-        else :	stop
+        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
-
+; 4 bands fit : IRAS 100 and 857 down to 353 GHz
+;fit_sed,"~/tmp/out/SED.RES",['IRAS','HFI'],iband=[0,0,0,1,1,1,1,1,1,0],beta=fixbeta,/noplot
+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
+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 
+endif
+
+print
+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)
+
+;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:*]
  
-;=== 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
+    print,'dustem_mpfit_run_vg: Current parameters values:',p_dim
   ENDIF
 ENDIF
  
@@ -100,7 +100,7 @@ endif else if keyword_set(albedo) then begin
 	ytit=textoidl('albedo')
 	plot,xrange,/nodata,xtit=xtit,ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,ys=ys,xr=xr,xsty=xsty
 	oplot,xrange,sca_tot/(sca_tot+abs_tot)
-	stop
+	
  
 endif else begin
  
@@ -116,7 +116,7 @@ endif else begin
  
 	Ngrains=n_tags(st.sed)-2
 	colors=dustem_grains_colors(Ngrains,ls,ps=ps)
-	
+;	stop
 	ytit=textoidl('\sigma_{ext}/N_H (cm^2/H)') 
 	if not keyword_set(ps) then tit='Dustem extinction'
  
@@ -18,7 +18,8 @@ endif else begin
 endelse
 endif
  
-st_grains=dustem_read_grain(dir_in+'data/'+(*!dustem_inputs).grain,/silent)
+;st_grains=dustem_read_grain(dir_in+'data/'+(*!dustem_inputs).grain,/silent)
+st_grains=dustem_read_grain(dir_in+'data/'+'GRAIN.DAT',/silent)
  
 ngrains = st_grains.ngrains
 colors=dustem_grains_colors(Ngrains,ls,ps=ps)
@@ -57,7 +57,7 @@ FOR i=1L,Ngrains DO BEGIN
   oplot,st.sed.wav,st.sed.(i)*fact,color=colors(i-1),psym=psym
 ENDFOR
 oplot,st.sed.wav,st.sed.(i)*fact,color=255,psym=psym
-stop
+
 defsysv,'!dustem_data',exists=exists
 if exists then begin
  
@@ -82,7 +82,8 @@ xtit=textoidl(&#39;\lambda (\mum)&#39;)
 ; Are observational data defined ?
 defsysv,'!dustem_data',exists=data_exists
  
-if keyword_set(print_ratio) then begin
+if ptr_valid(!dustem_data.polext) then begin
+   if keyword_set(print_ratio) then begin
  
 	; Correction de couleur
 	cc = 1.11
@@ -112,9 +113,14 @@ if keyword_set(print_ratio) then begin
 	tauV_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Vband)
         tauB = interpol(st.ext.ext_tot,st.ext.wav,Bband)
 	tauB_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Bband)
+    	print
+
  
 	; I353 and P353 are in fact nu_Pnu, in unit of erg/s/H, not erg/s/sr/H
-	; 1 W/m2/sr/1d20H = 1e-17 erg/s/sr/H
+                                ; 1 W/m2/sr/1d20H = 1e-17 erg/s/sr/H
+endif
+if ptr_valid(!dustem_data.polsed) then begin
+
         dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)
 	x=min(((*!dustem_data.polsed).wav-850)^2,jmin)
 	P353 = dustem_polsed(jmin)
@@ -132,26 +138,33 @@ if keyword_set(print_ratio) then begin
 	print
 	print,"Observational  constraint"
 	print,"-------------------------"
+	print,"I353/NH data        ",I353_data," MJy/sr for NH=1e20"
+	print,"P353/NH data        ",P353_data," MJy/sr for NH=1e20"
+        print,"P353/I353 data      ",P353_data/I353_data*100," %"
+        print
+        print,"Model results"
+	print,"-------------"
+        print,"I353/NH model       ",I353," MJy/sr for NH=1e20"
+	print,"P353/NH model       ",P353," MJy/sr for NH=1e20"
+	print,"P353/I353 model     ",P353/I353*100," %"
+
+        if keyword_set(print_ratio) then begin
+         print
+	print,"Observational  constraint"
+	print,"-------------------------"
+
         print,"tau_V data      ",tauV_data
         print,"(p/tau)_V data      ",pV_data/tauV_data*100," %"
         print,"p_V/A_V data        ",pV_data/(tauV_data*1.086)*100," %"
 	print,"with p_V data       ",pV_data," %"
-	print,"I353/NH data        ",I353_data," MJy/sr for NH=1e20"
-	print,"P353/NH data        ",P353_data," MJy/sr for NH=1e20"
-	print,"P353/I353 data      ",P353_data/I353_data*100," %"
         print,"Rs/v data           ",(P353_data/I353_data)/(pv_data/tauv_data)
         print,"RP/p data           ",(P353_data/1e20) / (pv_data*1e-21)," MJy/sr"
         print,"I353/Av data        ",(I353_data/1e20) / (1.086*tauv_data*1e-21)," MJy/sr"
 	print
-	print,"Model results"
-	print,"-------------"
         print,"(p/tau)_V model     ",pV/tauV*100," %"
         print,"p_V/A_V model       ",pV/(tauV*1.086)*100," %"
         print,"p_V/E(B-V) model       ",pV/((tauB-tauV)*1.086)*100," %"
 	print,"with p_V model      ",pV," %"
-	print,"I353/NH model       ",I353," MJy/sr for NH=1e20"
-	print,"P353/NH model       ",P353," MJy/sr for NH=1e20"
-	print,"P353/I353 model     ",P353/I353*100," %"
         print,"Rs/v model          ",P353/I353/(pV/tauV)
         print,"Rs/EBV model          ",P353/I353/(pV/(tauB-tauV))
         print,"RP/p model          ",(P353/1e20) / (pv*1e-21)," MJy/sr"
@@ -172,13 +185,14 @@ if keyword_set(print_ratio) then begin
         ;print,"Rs/v P353model/I353data/(pVmodel/tauVdata)	",P353/I353_data/(pV/tauV_data)
         ;print,"RP/p P353model/pVdata	",(P353/1e20)/(pV_data*1e-21)," MJy/sr"
         ;print
-
+     endif
+    endif    
  
 endif else if keyword_set(UV) then begin
  
 	if not keyword_set(ps) then tit='Dustem polarization cross-section'
 	ytit=textoidl('\sigma_{pol}/N_H (cm^2/H)')  
-
+if ptr_valid(!dustem_data.polext) then begin
 	if !dustem_show_plot EQ 2 or multi eq 2 then begin
 		;normpol = polext_tot 
 		normpol = st.polext.ext_tot 
@@ -193,7 +207,9 @@ endif else if keyword_set(UV) then begin
  
 	if multi eq 0 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog
 	if multi eq 1 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit='',ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)' 
-	if multi eq 2 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit='Normalized',tit='',xr=xr,yr=[-1,2],ylog=ylog,position=[0.17,0.14,0.95,0.35],/noerase,yticks=3,ymino=5,xticklen=0.1
+        if multi eq 2 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit='Normalized',tit='',xr=xr,yr=[-1,2],ylog=ylog,position=[0.17,0.14,0.95,0.35],/noerase,yticks=3,ymino=5,xticklen=0.1
+    
+
 	FOR i=0L,Ngrains-1 DO if aligned(i) then oplot,st.polext.wav,(st.polext.abs_grain(i)+st.polext.sca_grain(i))/normpol,color=colors(i+1),psym=psym,line=ls(i+1)
  
 	if total(aligned gt 0) then oplot,st.polext.wav,st.polext.ext_tot/normpol,psym=psym
@@ -203,8 +219,8 @@ endif else if keyword_set(UV) then begin
  		oplot,  (*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,ps=4
  		if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,yrms=(*!dustem_data.polext).sigma/dustem_polext
 	endif
-
-;	; Fit in (pmax,lmax,K)
+     endif
+                                ;	; Fit in (pmax,lmax,K)
 ;	x=st.polext.wav
 ;	logx=alog(x)
 ;
@@ -238,7 +254,7 @@ endif else if keyword_set(UV) then begin
 ;	; overplot the fit
 ;	;oplot,x(j),exp(res(0)+logx(j)*res(1)-K*logx(j)^2),col=150
  
-endif else if keyword_set(SED) then begin
+endif else if  ptr_valid(!dustem_data.polsed) then begin
  
         ;device, filename = 'plot008.eps', /color, /encapsulated
         if not keyword_set(ps) then tit='Polarization intensity in emission'
@@ -52,12 +52,13 @@ xtit=textoidl(&#39;a (\mu m)&#39;)
 norm_size = 1.d7
 norm_size = 1.d4
 xr /=1d3
+
 ;ytit='Mass log( 10!u29!n n!dH!u-1!n a dm/da ) '
 ytit=textoidl('4\pi/3 a^3dn/d ln a (cm^3/H)')
 if multi eq 0 then plot_oo, xr, /nodata,xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=''
 if multi eq 1 then plot_oo, xr, /nodata,xr=xr,/xs,xtit='',/ys,yr=yr,ytit=ytit,tit='',position=[0.17,0.45,0.95,0.95],xtickformat='(A1)' 
 colors=dustem_grains_colors(ntype,ls,ps=ps)
-
+;stop
 FOR i=0,ntype-1 DO begin
         ;if i le 1 and keyword_set(mergePAH) then begin
 	;	if i eq 0 then oplot, ax(0,0:nsize(0)-1)*1.e7, alog10(ava(0,0:nsize(0)-1)*yscl)+alog10(ava(1,0:nsize(1)-1)*yscl), line=ls(1),color = colors(1)
@@ -42,20 +42,22 @@ one_st={aligned:0.,law:&#39;&#39;,athresh:0.,plev:0.,pstiff:0.,TdsTg:0.,a0:0.,rvcut:0.}
 st=replicate(one_st,st_grains.Ngrains)
  
 st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool) 
-
+;count=0
 FOR i=0L,st_grains.Ngrains-1 DO BEGIN
-
-      IF st(i).aligned then begin
-
-	readf,unit,str,format='(A100)'
+   IF st(i).aligned then begin
+        ; if count EQ 0. then begin   
+         readf,unit,str,format='(A100)'   
+      ;endif 
+     
       	strv=str_sep(strcompress(strtrim(str,2)),' ')
       	Nstrv=n_elements(strv)
-	
+
 	ii=0L
-      	st(i).law=strv(ii) & ii=ii+1
+      	st(i).law=strv(ii) & ii=ii+1 
  
 	; IDG
-      	IF stregex(st(i).law, 'idg', /bool) THEN begin
+        IF stregex(st(i).law, 'idg', /bool) THEN begin
+           ii=ii+1 
       		st(i).TdsTg = strv(ii) & ii=ii+1
       		st(i).a0    = strv(ii) & ii=ii+1
       		st(i).rvcut = strv(ii) & ii=ii+1
@@ -66,19 +68,25 @@ FOR i=0L,st_grains.Ngrains-1 DO BEGIN
 	endif
  
 	; PARAMETRIC
-      	IF stregex(st(i).law, 'par', /bool) THEN begin
-		; Grain radius Threshold for alignment, given in microns and converted in cm
+        IF stregex(st(i).law, 'par', /bool) THEN begin
+                                ; Grain radius Threshold for
+                                ; alignment, given in microns and
+                                ; converted in cm
       		st(i).athresh = strv(ii) & ii=ii+1
       		st(i).pstiff  = strv(ii) & ii=ii+1
-      		st(i).plev    = strv(ii) & ii=ii+1
+                st(i).plev    = strv(ii) & ii=ii+1
+                ; count=1
 	endif
  
         IF !run_univ eq 1 then begin
 		st = replicate(st(i),st_grains.Ngrains)
-		st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool) 
+                st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool)
+                ;stop
 		break
 	endif
-
+       
+    
+       
       ENDIF
  
 ENDFOR
@@ -87,6 +95,7 @@ close,unit
 free_lun,unit
  
 full_st={keywords:key_str,anisG0:anisG0,gamma:0,grains:st}
+
 RETURN,full_st
  
 END
@@ -73,6 +73,7 @@ CASE !dustem_which OF
     file=dir_out+'out/'+'EXT.RES'
     st_ext=dustem_read_ext_lv(file,silent=silent)
 ;VGb : add polarization
+    
     if !run_pol then begin
         file=dir_out+'out/'+'SED_POL.RES'
         st_polsed=dustem_read_dustem_lv(file,silent=silent)
@@ -52,6 +52,7 @@ ENDIF
 ;=== clean the data structure
 !dustem_data.sed = ptr_new()
 !dustem_data.ext = ptr_new()
+
 if keyword_set(polext) then !dustem_data.polext = ptr_new()
 if keyword_set(polsed) then !dustem_data.polsed = ptr_new()
 if keyword_set(polfrac) then !dustem_data.polfrac = ptr_new()
@@ -108,7 +109,8 @@ IF keyword_set(polfrac) THEN BEGIN
 ENDIF
  
 IF keyword_set(polext) THEN BEGIN
-  wavs=polext.wave
+   wavs=polext.wave
+  
   ;=== Impose central wavelengths for photometric channels
   ind=where(polext.filter NE 'SPECTRUM',count)
   IF count NE 0 THEN BEGIN
@@ -116,7 +118,7 @@ IF keyword_set(polext) THEN BEGIN
   ENDIF
   ;=== define observations structure
   obs_st={instru_names:polext.instru,filt_names:polext.filter,wav:wavs, values:polext.spec,sigma:polext.error} 
-  ;=== fill !dustem_data
+                                ;=== fill !dustem_data
   !dustem_data.polext = ptr_new(obs_st)
 ENDIF
  
@@ -118,7 +118,7 @@
 ;  caller gracefully.  However, during the debugging process, it is
 ;  often useful to halt execution where the error occurred.  When you
 ;  set the NOCATCH keyword, MPFIT will not do any special error
-;  trapping, and execution will stop where ever the error occurred.
+;  trapping, and execution will stop whereever the error occurred.
 ;
 ;  MPFIT does not explicitly change the !ERROR_STATE variable
 ;  (although it may be changed implicitly if MPFIT calls MESSAGE).  It
@@ -164,9 +164,9 @@
 ;  individual values of PARINFO.MPSIDE.  Setting AUTODERIVATIVE=0 is
 ;  equivalent to resetting PARINFO.MPSIDE=3 for all parameters.
 ;
-;  However, be aware that even if the user requests explicit
-;  derivatives for some or all parameters, MPFIT will not always
-;  request explicit derivatives on ever user function call.
+;  Even if the user requests explicit derivatives for some or all
+;  parameters, MPFIT will not always request explicit derivatives on
+;  every user function call.
 ;
 ; EXPLICIT DERIVATIVES - CALLING INTERFACE
 ;
@@ -303,6 +303,14 @@
 ;  that the derivatives computed in two different ways have the same
 ;  numerical sign and the same order of magnitude, since these are the
 ;  most common programming mistakes.
+;
+;  A line of this form may also appear 
+;
+;   # FJAC_MASK = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
+;
+;  This line indicates for which parameters explicit derivatives are
+;  expected.  A list of all-1s indicates all explicit derivatives for
+;  all parameters are requested from the user function.
 ;    
 ;  
 ; CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD
@@ -448,6 +456,25 @@
 ;  constrained to be above 50.
 ;
 ;
+; RECURSION
+;
+;  Generally, recursion is not allowed.  As of version 1.77, MPFIT has
+;  recursion protection which does not allow a model function to
+;  itself call MPFIT.  Users who wish to perform multi-level
+;  optimization should investigate the 'EXTERNAL' function evaluation
+;  methods described below for hard-to-evaluate functions.  That
+;  method places more control in the user's hands.  The user can
+;  design a "recursive" application by taking care.
+;
+;  In most cases the recursion protection should be well-behaved.
+;  However, if the user is doing debugging, it is possible for the
+;  protection system to get "stuck."  In order to reset it, run the
+;  procedure:
+;     MPFIT_RESET_RECURSION
+;  and the protection system should get "unstuck."  It is save to call
+;  this procedure at any time.
+;
+;
 ; COMPATIBILITY
 ;
 ;  This function is designed to work with IDL 5.0 or greater.
@@ -457,6 +484,53 @@
 ;  of the IDL Virtual Machine.
 ;
 ;
+; DETERMINING THE VERSION OF MPFIT
+;
+;  MPFIT is a changing library.  Users of MPFIT may also depend on a
+;  specific version of the library being present.  As of version 1.70
+;  of MPFIT, a VERSION keyword has been added which allows the user to
+;  query which version is present.  The keyword works like this:
+;
+;    RESULT = MPFIT(/query, VERSION=version)
+;
+;  This call uses the /QUERY keyword to query the version number
+;  without performing any computations.  Users of MPFIT can call this
+;  method to determine which version is in the IDL path before
+;  actually using MPFIT to do any numerical work.  Upon return, the
+;  VERSION keyword contains the version number of MPFIT, expressed as
+;  a string of the form 'X.Y' where X and Y are integers.
+;
+;  Users can perform their own version checking, or use the built-in
+;  error checking of MPFIT.  The MIN_VERSION keyword enforces the
+;  requested minimum version number.  For example,
+;
+;    RESULT = MPFIT(/query, VERSION=version, MIN_VERSION='1.70')
+;
+;  will check whether the accessed version is 1.70 or greater, without
+;  performing any numerical processing.
+;
+;  The VERSION and MIN_VERSION keywords were added in MPFIT
+;  version 1.70 and later.  If the caller attempts to use the VERSION
+;  or MIN_VERSION keywords, and an *older* version of the code is
+;  present in the caller's path, then IDL will throw an 'unknown
+;  keyword' error.  Therefore, in order to be robust, the caller, must
+;  use exception handling.  Here is an example demanding at least
+;  version 1.70.
+;
+;    MPFIT_OK = 0  & VERSION = '<unknown>'
+;    CATCH, CATCHERR
+;    IF CATCHERR EQ 0 THEN MPFIT_OK = MPFIT(/query, VERSION=version, $
+;                                         MIN_VERSION='1.70')
+;    CATCH, /CANCEL
+;
+;    IF NOT MPFIT_OK THEN $
+;      MESSAGE, 'ERROR: you must have MPFIT version 1.70 or higher in '+$
+;             'your path (found version '+version+')'
+;
+;  Of course, the caller can also do its own version number
+;  requirements checking.
+;
+;
 ; HARD-TO-COMPUTE FUNCTIONS: "EXTERNAL" EVALUATION
 ;
 ;  The normal mode of operation for MPFIT is for the user to pass a
@@ -525,6 +599,8 @@
 ; RETURNS:
 ;
 ;   Returns the array of best-fit parameters.
+;   Exceptions: 
+;      * if /QUERY is set (see QUERY).
 ;
 ;
 ; KEYWORD PARAMETERS:
@@ -537,8 +613,31 @@
 ;                    NOTE: to supply your own explicit derivatives,
 ;                      explicitly pass AUTODERIVATIVE=0
 ;
-;   BESTNORM - the value of the summed squared weighted residuals for
-;              the returned parameter values, i.e. TOTAL(DEVIATES^2).
+;   BESTNORM - upon return, the value of the summed squared weighted
+;              residuals for the returned parameter values,
+;              i.e. TOTAL(DEVIATES^2).
+;
+;   BEST_FJAC - upon return, BEST_FJAC contains the Jacobian, or
+;               partial derivative, matrix for the best-fit model.
+;               The values are an array,
+;               ARRAY(N_ELEMENTS(DEVIATES),NFREE) where NFREE is the
+;               number of free parameters.  This array is only
+;               computed if /CALC_FJAC is set, otherwise BEST_FJAC is
+;               undefined.
+;
+;               The returned array is such that BEST_FJAC[I,J] is the
+;               partial derivative of DEVIATES[I] with respect to
+;               parameter PARMS[PFREE_INDEX[J]].  Note that since
+;               deviates are (data-model)*weight, the Jacobian of the
+;               *deviates* will have the opposite sign from the
+;               Jacobian of the *model*, and may be scaled by a
+;               factor.
+;
+;   BEST_RESID - upon return, an array of best-fit deviates.
+;
+;   CALC_FJAC - if set, then calculate the Jacobian and return it in
+;               BEST_FJAC.  If not set, then the return value of
+;               BEST_FJAC is undefined.
 ;
 ;   COVAR - the covariance matrix for the set of parameters returned
 ;           by MPFIT.  The matrix is NxN where N is the number of
@@ -548,9 +647,11 @@
 ;           Parameter errors are also returned in PERROR.
 ;
 ;           To compute the correlation matrix, PCOR, use this example:
-;           IDL> PCOR = COV * 0
-;           IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $
-;                PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j])
+;                  PCOR = COV * 0
+;                  FOR i = 0, n-1 DO FOR j = 0, n-1 DO $
+;                    PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j])
+;           or equivalently, in vector notation,
+;                  PCOR = COV / (PERROR # PERROR)
 ;
 ;           If NOCOVAR is set or MPFIT terminated abnormally, then
 ;           COVAR is set to a scalar with value !VALUES.D_NAN.
@@ -558,7 +659,8 @@
 ;   DOF - number of degrees of freedom, computed as
 ;             DOF = N_ELEMENTS(DEVIATES) - NFREE
 ;         Note that this doesn't account for pegged parameters (see
-;         NPEGGED).
+;         NPEGGED).  It also does not account for data points which
+;         are assigned zero weight by the user function.
 ;
 ;   ERRMSG - a string error or warning message is returned.
 ;
@@ -699,10 +801,24 @@
 ;                 change to a printable character like 'q'.
 ;
 ;   MAXITER - The maximum number of iterations to perform.  If the
-;             number is exceeded, then the STATUS value is set to 5
-;             and MPFIT returns.
+;             number of calculation iterations exceeds MAXITER, then
+;             the STATUS value is set to 5 and MPFIT returns.  
+;
+;             If MAXITER EQ 0, then MPFIT does not iterate to adjust
+;             parameter values; however, the user function is evaluated
+;             and parameter errors/covariance/Jacobian are estimated
+;             before returning.
 ;             Default: 200 iterations
 ;
+;   MIN_VERSION - The minimum requested version number.  This must be
+;                 a scalar string of the form returned by the VERSION
+;                 keyword.  If the current version of MPFIT does not
+;                 satisfy the minimum requested version number, then,
+;                    MPFIT(/query, min_version='...') returns 0
+;                    MPFIT(...) returns NAN
+;                 Default: no version number check
+;                 NOTE: MIN_VERSION was added in MPFIT version 1.70
+;
 ;   NFEV - the number of MYFUNCT function evaluations performed.
 ;
 ;   NFREE - the number of free parameters in the fit.  This includes
@@ -761,6 +877,19 @@
 ;              DOF     = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom
 ;              PCERROR = PERROR * SQRT(BESTNORM / DOF)   ; scaled uncertainties
 ;
+;   PFREE_INDEX - upon return, PFREE_INDEX contains an index array
+;                 which indicates which parameter were allowed to
+;                 vary.  I.e. of all the parameters PARMS, only
+;                 PARMS[PFREE_INDEX] were varied.
+;
+;   QUERY - if set, then MPFIT() will return immediately with one of
+;           the following values:
+;                 1 - if MIN_VERSION is not set
+;                 1 - if MIN_VERSION is set and MPFIT satisfies the minimum
+;                 0 - if MIN_VERSION is set and MPFIT does not satisfy it
+;           The VERSION output keyword is always set upon return.
+;           Default: not set.
+;
 ;   QUIET - set this keyword when no textual output should be printed
 ;           by MPFIT
 ;
@@ -827,6 +956,14 @@
 ;             the function and its derivatives.  This status indicator
 ;             is neither an error nor a convergence indicator.
 ;
+;   VERSION - upon return, VERSION will be set to the MPFIT internal
+;             version number.  The version number will be a string of
+;             the form "X.Y" where X is a major revision number and Y
+;             is a minor revision number.
+;             NOTE: the VERSION keyword was not present before 
+;               MPFIT version number 1.70, therefore, callers must 
+;               use exception handling when using this keyword.
+;
 ;   XTOL - a nonnegative input variable. Termination occurs when the
 ;          relative error between two consecutive iterates is at most
 ;          XTOL (and STATUS is accordingly set to 2 or 3).  Therefore,
@@ -965,18 +1102,24 @@
 ;     with MPFIT," in proc. Astronomical Data Analysis Software and
 ;     Systems XVIII, Quebec, Canada, ASP Conference Series, Vol. XXX, eds.
 ;     D. Bohlender, P. Dowler & D. Durand (Astronomical Society of the
-;     Pacific: San Francisco), p. XXX
-;       http://arxiv.org/XXXXX
+;     Pacific: San Francisco), p. 251-254 (ISBN: 978-1-58381-702-5)
+;       http://arxiv.org/abs/0902.2850
+;       Link to NASA ADS: http://adsabs.harvard.edu/abs/2009ASPC..411..251M
+;       Link to ASP: http://aspbooks.org/a/volumes/table_of_contents/411
+;
+;   Refer to the MPFIT website as:
+;       http://purl.com/net/mpfit
 ;
-;   MINPACK-1, Jorge More', available from netlib.
+;   MINPACK-1 software, by Jorge More' et al, available from netlib.
 ;     http://www.netlib.org/
 ;
 ;   "Optimization Software Guide," Jorge More' and Stephen Wright, 
 ;     SIAM, *Frontiers in Applied Mathematics*, Number 14.
+;     (ISBN: 978-0-898713-22-0)
 ;
-;   More', J. 1977, ``The Levenberg-Marquardt Algorithm: 
-;     Implementation and Theory,'' in *Numerical Analysis*, 
-;     vol. 630, ed. G. A. Watson (Springer-Verlag: Berlin), 105
+;   More', J. 1978, "The Levenberg-Marquardt Algorithm: Implementation
+;     and Theory," in Numerical Analysis, vol. 630, ed. G. A. Watson
+;     (Springer-Verlag: Berlin), p. 105 (DOI: 10.1007/BFb0067690 )
 ;
 ; MODIFICATION HISTORY:
 ;   Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM
@@ -1141,12 +1284,33 @@
 ;      the covariance matrix to be computed, CM, 14 Apr 2009
 ;   Avoid numerical underflow while solving for the LM parameter,
 ;      (thanks to Sergey Koposov) CM, 14 Apr 2009
-;
-;  $Id: mpfit.pro,v 1.64 2009/05/05 05:37:49 craigm Exp $
+;   Use individual functions for all possible MPFIT_CALL permutations,
+;      (and make sure the syntax is right) CM, 01 Sep 2009
+;   Correct behavior of MPMAXSTEP when some parameters are frozen,
+;      thanks to Josh Destree, CM, 22 Nov 2009
+;   Update the references section, CM, 22 Nov 2009
+;   1.70 - Add the VERSION and MIN_VERSION keywords, CM, 22 Nov 2009
+;   1.71 - Store pre-calculated revision in common, CM, 23 Nov 2009
+;   1.72-1.74 - Documented alternate method to compute correlation matrix,
+;          CM, 05 Feb 2010
+;   1.75 - Enforce TIED constraints when preparing to terminate the
+;          routine, CM, 2010-06-22
+;   1.76 - Documented input keywords now are not modified upon output,
+;          CM, 2010-07-13
+;   1.77 - Upon user request (/CALC_FJAC), compute Jacobian matrix and
+;          return in BEST_FJAC; also return best residuals in
+;          BEST_RESID; also return an index list of free parameters as
+;          PFREE_INDEX; add a fencepost to prevent recursion
+;          CM, 2010-10-27
+;   1.79 - Documentation corrections.  CM, 2011-08-26
+;   1.81 - Fix bug in interaction of AUTODERIVATIVE=0 and .MPSIDE=3;
+;          Document FJAC_MASK. CM, 2012-05-08
+;
+;  $Id: mpfit.pro,v 1.82 2012/09/27 23:59:44 cmarkwar Exp $
 ;-
 ; Original MINPACK by More' Garbow and Hillstrom, translated with permission
 ; Modifications and enhancements are:
-; Copyright (C) 1997-2003, 2004, 2005, 2006, 2007, 2008, 2009, Craig Markwardt
+; Copyright (C) 1997-2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Craig Markwardt
 ; This software is provided as is without any warranty whatsoever.
 ; Permission to use, copy, modify, and distribute modified or
 ; unmodified copies is granted, provided this copyright and disclaimer
@@ -1215,6 +1379,44 @@ pro mpfit_setmachar, double=isdouble
 end
  
  
+; Call user function with no _EXTRA parameters
+function mpfit_call_func_noextra, fcn, x, fjac, _EXTRA=extra
+  if n_params() EQ 2 then begin
+     return, call_function(fcn, x)
+  endif else begin
+     return, call_function(fcn, x, fjac)
+  endelse
+end
+
+; Call user function with _EXTRA parameters
+function mpfit_call_func_extra, fcn, x, fjac, _EXTRA=extra
+  if n_params() EQ 2 then begin
+     return, call_function(fcn, x, _EXTRA=extra)
+  endif else begin
+     return, call_function(fcn, x, fjac, _EXTRA=extra)
+  endelse
+end
+
+; Call user procedure with no _EXTRA parameters
+function mpfit_call_pro_noextra, fcn, x, fjac, _EXTRA=extra
+  if n_params() EQ 2 then begin
+     call_procedure, fcn, x, f
+  endif else begin
+     call_procedure, fcn, x, f, fjac
+  endelse
+  return, f
+end
+
+; Call user procedure with _EXTRA parameters
+function mpfit_call_pro_extra, fcn, x, fjac, _EXTRA=extra
+  if n_params() EQ 2 then begin
+     call_procedure, fcn, x, f, _EXTRA=extra
+  endif else begin
+     call_procedure, fcn, x, f, fjac, _EXTRA=extra
+  endelse
+  return, f
+end
+
  
 ;; Call user function or procedure, with _EXTRA or not, with
 ;; derivatives or not.
@@ -1224,36 +1426,19 @@ function mpfit_call, fcn, x, fjac, _EXTRA=extra
   common mpfit_config, mpconfig
  
   if keyword_set(mpconfig.qanytied) then mpfit_tie, x, mpconfig.ptied
-    
-  ;; Decide whether we are calling a procedure or function
-  if mpconfig.proc then proc = 1 else proc = 0
-  mpconfig.nfev = mpconfig.nfev + 1
  
-  if proc then begin
-      if n_params() EQ 3 then begin
-          if n_elements(extra) GT 0 then $
-            call_procedure, fcn, x, f, fjac, _EXTRA=extra $
-          else $
-            call_procedure, fcn, x, f, fjac
-      endif else begin
-          if n_elements(extra) GT 0 then $
-            call_procedure, fcn, x, f, _EXTRA=extra $
-          else $
-            call_procedure, fcn, x, f
-      endelse
+  ;; Decide whether we are calling a procedure or function, and 
+  ;; with/without FUNCTARGS
+  proname = 'MPFIT_CALL'
+  proname = proname + ((mpconfig.proc) ? '_PRO' : '_FUNC')
+  proname = proname + ((n_elements(extra) GT 0) ? '_EXTRA' : '_NOEXTRA')
+
+  if n_params() EQ 2 then begin
+     f = call_function(proname, fcn, x, _EXTRA=extra)
   endif else begin
-      if n_params() EQ 3 then begin
-          if n_elements(extra) GT 0 then $
-            f = call_function(fcn, x, fjac, _EXTRA=extra) $
-          else $
-            f = call_function(fcn, x, fjac)
-      endif else begin
-          if n_elements(extra) GT 0 then $
-            f = call_function(fcn, x, _EXTRA=extra) $
-          else $
-            f = call_function(fcn, x)
-      endelse
-  endelse  
+     f = call_function(proname, fcn, x, fjac, _EXTRA=extra)
+  endelse
+  mpconfig.nfev = mpconfig.nfev + 1
  
   if n_params() EQ 2 AND mpconfig.damp GT 0 then begin
       damp = mpconfig.damp[0]
@@ -1309,12 +1494,14 @@ function mpfit_fdjac2, fcn, x, fvec, step, ulimited, ulimit, dside, $
   ;;   2. AUTODERIVATIVE=1, but P[i].MPSIDE EQ 3
  
   if keyword_set(autoderiv) EQ 0 OR max(dside[ifree] EQ 3) EQ 1 then begin
-      fjac = intarr(nall)
+      fjac_mask = intarr(nall)
       ;; Specify which parameters need derivatives
-      ;;            ---- Case 2 ------     ----- Case 1 -----
-      fjac[ifree] = (dside[ifree] EQ 3) OR (keyword_set(autoderiv) EQ 0)
-      if has_debug_deriv then print, fjac, format='("# FJAC_MASK = ",100000(I0," ",:))'
+      ;;                  ---- Case 2 ------     ----- Case 1 -----
+      fjac_mask[ifree] = (dside[ifree] EQ 3) OR (keyword_set(autoderiv) EQ 0)
+      if has_debug_deriv then $
+        print, fjac_mask, format='("# FJAC_MASK = ",100000(I0," ",:))'
  
+      fjac = fjac_mask  ;; Pass the mask to the calling function as FJAC
       mperr = 0
       fp = mpfit_call(fcn, xall, fjac, _EXTRA=fcnargs)
       iflag = mperr
@@ -1340,12 +1527,13 @@ function mpfit_fdjac2, fcn, x, fvec, step, ulimited, ulimit, dside, $
       ;; Select only the free parameters
       if n_elements(ifree) LT nall then $
         fjac = reform(fjac[*,ifree], m, n, /overwrite)
-
-      ;; If there are no more free parameters to analyze, then
-      ;; return now, (but not if we are debugging the derivatives)
-      if ((keyword_set(autoderiv) EQ 0) OR $
-          (min(dside[ifree]) EQ 1) OR $
-          (has_debug_deriv EQ 0)) then return, fjac
+      
+      ;; Are we done computing derivatives?  The answer is, YES, if we
+      ;; computed explicit derivatives for all free parameters, EXCEPT
+      ;; when we are going on to compute debugging derivatives.
+      if min(fjac_mask[ifree]) EQ 1 AND NOT has_debug_deriv then begin
+          return, fjac
+      endif
   endif
  
   ;; Final output array, if it was not already created above
@@ -2174,7 +2362,7 @@ pro mpfit_defiter, fcn, x, iter, fnorm, FUNCTARGS=fcnargs, $
   if n_elements(fmt) GT 0 then begin
       call_procedure, iterprint, x, format=fmt, _EXTRA=iterargs
   endif else begin
-      if n_elements(pformat) EQ 0 then pformat = '(G20.6)'
+      if n_elements(pformat) EQ 0 then pformat = '(G40.6)'
       parname = 'P('+strtrim(iprint,2)+')'
       pformats = strarr(nprint) + pformat
  
@@ -2398,6 +2586,45 @@ function mpfit_covar, rr, ipvt, tol=tol
   return, r
 end
  
+;; Parse the RCSID revision number
+function mpfit_revision
+  ;; NOTE: this string is changed every time an RCS check-in occurs
+  revision = '$Revision: 1.82 $'
+
+  ;; Parse just the version number portion
+  revision = stregex(revision,'\$'+'Revision: *([0-9.]+) *'+'\$',/extract,/sub)
+  revision = revision[1]
+  return, revision
+end
+
+;; Parse version numbers of the form 'X.Y'
+function mpfit_parse_version, version
+  sz = size(version)
+  if sz[sz[0]+1] NE 7 then return, 0
+
+  s = stregex(version[0], '^([0-9]+)\.([0-9]+)$', /extract,/sub) 
+  if s[0] NE version[0] then return, 0
+  return, long(s[1:2])
+end
+
+;; Enforce a minimum version number
+function mpfit_min_version, version, min_version
+  mv = mpfit_parse_version(min_version)
+  if mv[0] EQ 0 then return, 0
+  v  = mpfit_parse_version(version)
+
+  ;; Compare version components
+  if v[0] LT mv[0] then return, 0
+  if v[1] LT mv[1] then return, 0
+  return, 1
+end
+
+; Manually reset recursion fencepost if the user gets in trouble
+pro mpfit_reset_recursion
+  common mpfit_fencepost, mpfit_fencepost_active
+  mpfit_fencepost_active = 0
+end
+
 ;     **********
 ;
 ;     subroutine lmdif
@@ -2579,24 +2806,48 @@ end
 ;
 ;     **********
 function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
-                ftol=ftol, xtol=xtol, gtol=gtol, epsfcn=epsfcn, resdamp=damp, $
+                ftol=ftol0, xtol=xtol0, gtol=gtol0, epsfcn=epsfcn, $
+                resdamp=damp0, $
                 nfev=nfev, maxiter=maxiter, errmsg=errmsg, $
-                factor=factor, nprint=nprint, STATUS=info, $
-                iterproc=iterproc, iterargs=iterargs, iterstop=ss,$
+                factor=factor0, nprint=nprint0, STATUS=info, $
+                iterproc=iterproc0, iterargs=iterargs, iterstop=ss,$
                 iterkeystop=iterkeystop, $
                 niter=iter, nfree=nfree, npegged=npegged, dof=dof, $
-                diag=diag, rescale=rescale, autoderivative=autoderiv, $
-                perror=perror, covar=covar, nocovar=nocovar, bestnorm=fnorm, $
+                diag=diag, rescale=rescale, autoderivative=autoderiv0, $
+                pfree_index=ifree, $
+                perror=perror, covar=covar, nocovar=nocovar, $
+                bestnorm=fnorm, best_resid=fvec, $
+                best_fjac=output_fjac, calc_fjac=calc_fjac, $
                 parinfo=parinfo, quiet=quiet, nocatch=nocatch, $
-                fastnorm=fastnorm, proc=proc, query=query, $
+                fastnorm=fastnorm0, proc=proc, query=query, $
                 external_state=state, external_init=extinit, $
-                external_fvec=efvec, external_fjac=efjac
-
+                external_fvec=efvec, external_fjac=efjac, $
+                version=version, min_version=min_version0
+  
+  
   COMPILE_OPT strictarr
   info = 0L
   errmsg = ''
  
-  if keyword_set(query) then return, 1
+  ;; Compute the revision number, to be returned in the VERSION and
+  ;; QUERY keywords.
+  common mpfit_revision_common, mpfit_revision_str
+  if n_elements(mpfit_revision_str) EQ 0 then $
+     mpfit_revision_str = mpfit_revision()
+  version = mpfit_revision_str
+
+  if keyword_set(query) then begin
+     if n_elements(min_version0) GT 0 then $
+        if mpfit_min_version(version, min_version0[0]) EQ 0 then $
+           return, 0
+     return, 1
+  endif
+
+  if n_elements(min_version0) GT 0 then $
+     if mpfit_min_version(version, min_version0[0]) EQ 0 then begin
+     message, 'ERROR: minimum required version '+min_version0[0]+' not satisfied', /info
+     return, !values.d_nan
+  endif
  
   if n_params() EQ 0 then begin
       message, "USAGE: PARMS = MPFIT('MYFUNCT', START_PARAMS, ... )", /info
@@ -2605,15 +2856,15 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
  
   ;; Use of double here not a problem since f/x/gtol are all only used
   ;; in comparisons
-  if n_elements(ftol) EQ 0 then ftol = 1.D-10
-  if n_elements(xtol) EQ 0 then xtol = 1.D-10
-  if n_elements(gtol) EQ 0 then gtol = 1.D-10
-  if n_elements(factor) EQ 0 then factor = 100.
-  if n_elements(nprint) EQ 0 then nprint = 1
-  if n_elements(iterproc) EQ 0 then iterproc = 'MPFIT_DEFITER'
-  if n_elements(autoderiv) EQ 0 then autoderiv = 1
-  if n_elements(fastnorm) EQ 0 then fastnorm = 0
-  if n_elements(damp) EQ 0 then damp = 0 else damp = damp[0]
+  if n_elements(ftol0) EQ 0 then ftol = 1.D-10 else ftol = ftol0[0]
+  if n_elements(xtol0) EQ 0 then xtol = 1.D-10 else xtol = xtol0[0]
+  if n_elements(gtol0) EQ 0 then gtol = 1.D-10 else gtol = gtol0[0]
+  if n_elements(factor0) EQ 0 then factor = 100. else factor = factor0[0]
+  if n_elements(nprint0) EQ 0 then nprint = 1 else nprint = nprint0[0]
+  if n_elements(iterproc0) EQ 0 then iterproc = 'MPFIT_DEFITER' else iterproc = iterproc0[0]
+  if n_elements(autoderiv0) EQ 0 then autoderiv = 1 else autoderiv = autoderiv0[0]
+  if n_elements(fastnorm0) EQ 0 then fastnorm = 0 else fastnorm = fastnorm0[0]
+  if n_elements(damp0) EQ 0 then damp = 0 else damp = damp0[0]
  
   ;; These are special configuration parameters that can't be easily
   ;; passed by MPFIT directly.
@@ -2631,6 +2882,18 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
   nfree = 0L
   npegged = 0L
   dof = 0L
+  output_fjac = 0L
+
+  ;; Set up a persistent fencepost that prevents recursive calls
+  common mpfit_fencepost, mpfit_fencepost_active
+  if n_elements(mpfit_fencepost_active) EQ 0 then mpfit_fencepost_active = 0
+  if mpfit_fencepost_active then begin
+      errmsg = 'ERROR: recursion detected; you cannot run MPFIT recursively'
+      goto, TERMINATE
+  endif
+  ;; Only activate the fencepost if we are not in debugging mode
+  if NOT keyword_set(nocatch) then mpfit_fencepost_active = 1
+
  
   ;; Parameter damping doesn't work when user is providing their own
   ;; gradients.
@@ -2669,8 +2932,9 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
   ;; Handle error conditions gracefully
   if NOT keyword_set(nocatch) then begin
       catch, catcherror
-      if catcherror NE 0 then begin
+      if catcherror NE 0 then begin  ;; An error occurred!!!
           catch, /cancel
+          mpfit_fencepost_active = 0
           err_string = ''+!error_state.msg
           message, /cont, 'Error detected while '+catch_msg+':'
           message, /cont,    err_string
@@ -2795,8 +3059,6 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
       errmsg = 'ERROR: MPMINSTEP is greater than MPMAXSTEP'
       goto, TERMINATE
   endif
-  wh = where(qmin OR qmax, ct)
-  qminmax = ct GT 0
  
   ;; Finish up the free parameters
   ifree = where(pfixed NE 1, nfree)
@@ -2819,6 +3081,12 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
   ;; Compose only VARYING parameters
   xnew = xall      ;; xnew is the set of parameters to be returned
   x = xnew[ifree]  ;; x is the set of free parameters
+  ; Same for min/max step diagnostics
+  qmin = qmin[ifree]  & minstep = minstep[ifree]
+  qmax = qmax[ifree]  & maxstep = maxstep[ifree]
+  wh = where(qmin OR qmax, ct)
+  qminmax = ct GT 0
+
  
   ;; LIMITED parameters ?
   mpfit_parinfo, parinfo, tagnames, 'LIMITED', limited, status=st1
@@ -3012,6 +3280,8 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
   iflag = 2
   if NOT isext then begin
       catch_msg = 'calling MPFIT_FDJAC2'
+      ;; NOTE!  If you change this call then change the one during
+      ;; clean-up as well!
       fjac = mpfit_fdjac2(fcn, x, fvec, step, qulim, ulim, dside, $
                           iflag=iflag, epsfcn=epsfcn, $
                           autoderiv=autoderiv, dstep=dstep, $
@@ -3043,6 +3313,7 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
       ;; See if any "pegged" values should keep their derivatives
       if (nlpeg GT 0) then begin
           ;; Total derivative of sum wrt lower pegged parameters
+          ;;   Note: total(fvec*fjac[*,i]) is d(CHI^2)/dX[i]
           for i = 0L, nlpeg-1 do begin
               sum = total(fvec * fjac[*,whlpeg[i]])
               if sum GT 0 then fjac[*,whlpeg[i]] = 0
@@ -3057,24 +3328,42 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
       endif
   endif
  
+  ;; Save a copy of the Jacobian if the user requests it...
+  if keyword_set(calc_fjac) then output_fjac = fjac
+
+  ;; =====================
   ;; Compute the QR factorization of the jacobian
   catch_msg = 'calling MPFIT_QRFAC'
-  mpfit_qrfac, fjac, ipvt, wa1, wa2, /pivot
-
+  ;;  IN:      Jacobian        
+  ;; OUT:      Hh Vects  Permutation  RDIAG  ACNORM
+  mpfit_qrfac, fjac,     ipvt,        wa1,   wa2, /pivot
+  ;; Jacobian - jacobian matrix computed by mpfit_fdjac2
+  ;; Hh vects - house holder vectors from QR factorization & R matrix
+  ;; Permutation - permutation vector for pivoting
+  ;; RDIAG - diagonal elements of R matrix
+  ;; ACNORM - norms of input Jacobian matrix before factoring
+
+  ;; =====================
   ;; On the first iteration if "diag" is unspecified, scale
   ;; according to the norms of the columns of the initial jacobian
   catch_msg = 'rescaling diagonal elements'
   if (iter EQ 1) then begin
+      ;; Input: WA2 = root sum of squares of original Jacobian matrix
+      ;;        DIAG = user-requested diagonal (not documented!)
+      ;;        FACTOR = user-requested norm factor (not documented!)
+      ;; Output: DIAG = Diagonal scaling values
+      ;;         XNORM = sum of squared scaled residuals
+      ;;         DELTA = rescaled XNORM
  
       if NOT keyword_set(rescale) OR (n_elements(diag) LT n) then begin
-          diag = wa2
-          wh = where (diag EQ 0, ct)
+          diag = wa2 ;; Calculated from original Jacobian
+          wh = where (diag EQ 0, ct) ;; Handle zero values
           if ct GT 0 then diag[wh] = one
       endif
  
       ;; On the first iteration, calculate the norm of the scaled x
       ;; and initialize the step bound delta 
-      wa3 = diag * x
+      wa3 = diag * x           ;; WA3 is temp variable
       xnorm = mpfit_enorm(wa3)
       delta = factor*xnorm
       if delta EQ zero then delta = zero + factor
@@ -3176,7 +3465,7 @@ function mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
           nwa1 = wa1 * alpha
           whmax = where(qmax AND maxstep GT 0, ct)
           if ct GT 0 then begin
-              mrat = max(abs(nwa1[whmax])/abs(maxstep[ifree[whmax]]))
+              mrat = max(abs(nwa1[whmax])/abs(maxstep[whmax]))
               if mrat GT 1 then alpha = alpha / mrat
           endif
       endif          
@@ -3316,6 +3605,7 @@ TERMINATE:
   iflag = 0
   if n_elements(xnew) EQ 0 then goto, FINAL_RETURN
   if nfree EQ 0 then xnew = xall else xnew[ifree] = x
+  if n_elements(qanytied) GT 0 then if qanytied then mpfit_tie, xnew, ptied
   dof = n_elements(fvec) - nfree
  
  
@@ -3350,7 +3640,8 @@ TERMINATE:
                  (qllim AND (x EQ llim)), npegged)
   endif
  
-  if fcn NE '(EXTERNAL)' AND nprint GT 0 AND info GT 0 then begin
+  ;; Calculate final function value (FNORM) and residuals (FVEC)
+  if isext EQ 0 AND nprint GT 0 AND info GT 0 then begin
       catch_msg = 'calling '+fcn
       fvec = mpfit_call(fcn, xnew, _EXTRA=fcnargs)
       catch_msg = 'in the termination phase'
@@ -3394,10 +3685,12 @@ TERMINATE:
             perror[wh] = sqrt(covar[wh, wh])
       endif
   endif
+
 ;  catch_msg = 'returning the result'
 ;  profvals.mpfit = profvals.mpfit + (systime(1) - prof_start)
  
   FINAL_RETURN:
+  mpfit_fencepost_active = 0
   nfev = mpconfig.nfev
   if n_elements(xnew) EQ 0 then return, !values.d_nan
   return, xnew
@@ -73,12 +73,61 @@
 ;  using the ERROR_CODE common block variable, as described
 ;  below under the MPFIT_ERROR common block definition.
 ;
-;  See the discussion under "EXPLICIT DERIVATIVES" and AUTODERIVATIVE
-;  in MPFIT.PRO if you wish to compute the derivatives for yourself.
-;  AUTODERIVATIVE is accepted by MPFITFUN and passed directly to
-;  MPFIT.  The user function must accept one additional parameter, DP,
-;  which contains the derivative of the user function with respect to
-;  each parameter at each data point, as described in MPFIT.PRO.
+;  MPFIT by default calculates derivatives numerically via a finite
+;  difference approximation.  However, the user function *may*
+;  calculate the derivatives if desired, but only if the model
+;  function is declared with an additional position parameter, DP, as
+;  described below.
+;
+;  To enable explicit derivatives for all parameters, set
+;  AUTODERIVATIVE=0.
+;
+;  When AUTODERIVATIVE=0, the user function is responsible for
+;  calculating the derivatives of the user function with respect to
+;  each parameter.  The user function should be declared as follows:
+;
+;    ;
+;    ; MYFUNCT - example user function
+;    ;   P - input parameter values (N-element array)
+;    ;   DP - upon input, an N-vector indicating which parameters
+;    ;          to compute derivatives for; 
+;    ;        upon output, the user function must return
+;    ;          an ARRAY(M,N) of derivatives in this keyword
+;    ;   (keywords) - any other keywords specified by FUNCTARGS
+;    ; RETURNS - function values
+;    ;
+;    FUNCTION MYFUNCT, x, p, dp [, (additional keywords if desired)]
+;     model = F(x, p)         ;; Model function
+;     
+;     if n_params() GT 2 then begin
+;       ; Create derivative and compute derivative array
+;       requested = dp   ; Save original value of DP
+;       dp = make_array(n_elements(x), n_elements(p), value=x[0]*0)
+;
+;       ; Compute derivative if requested by caller
+;       for i = 0, n_elements(p)-1 do if requested(i) NE 0 then $
+;         dp(*,i) = FGRAD(x, p, i)
+;     endif
+;    
+;     return, resid
+;    END
+;
+;  where FGRAD(x, p, i) is a model function which computes the
+;  derivative of the model F(x,p) with respect to parameter P(i) at X.
+;
+;  Derivatives should be returned in the DP array. DP should be an
+;  ARRAY(m,n) array, where m is the number of data points and n is the
+;  number of parameters.  DP[i,j] is the derivative of the ith
+;  function value with respect to the jth parameter.
+;
+;  MPFIT may not always request derivatives from the user function.
+;  In those cases, the parameter DP is not passed.  Therefore
+;  functions can use N_PARAMS() to indicate whether they must compute
+;  the derivatives or not.
+;
+;  For additional information about explicit derivatives, including
+;  additional settings and debugging options, see the discussion under
+;  "EXPLICIT DERIVATIVES" and AUTODERIVATIVE in MPFIT.PRO.
 ;
 ; CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD
 ;
@@ -253,6 +302,21 @@
 ;   BESTNORM - the value of the summed squared residuals for the
 ;              returned parameter values.
 ;
+;   BEST_FJAC - upon return, BEST_FJAC contains the Jacobian, or
+;               partial derivative, matrix for the best-fit model.
+;               The values are an array,
+;               ARRAY(N_ELEMENTS(DEVIATES),NFREE) where NFREE is the
+;               number of free parameters.  This array is only
+;               computed if /CALC_FJAC is set, otherwise BEST_FJAC is
+;               undefined.
+;
+;               The returned array is such that BEST_FJAC[I,J] is the
+;               partial derivative of the model with respect to
+;               parameter PARMS[PFREE_INDEX[J]].
+;
+;   BEST_RESID - upon return, an array of best-fit deviates,
+;                normalized by the weights or errors.
+;
 ;   COVAR - the covariance matrix for the set of parameters returned
 ;           by MPFIT.  The matrix is NxN where N is the number of
 ;           parameters.  The square root of the diagonal elements
@@ -260,10 +324,12 @@
 ;           parameters IF errors were treated "properly" in MYFUNC.
 ;           Parameter errors are also returned in PERROR.
 ;
-;           To compute the correlation matrix, PCOR, use this:
-;           IDL> PCOR = COV * 0
-;           IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $
-;                PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j])
+;           To compute the correlation matrix, PCOR, use this example:
+;                  PCOR = COV * 0
+;                  FOR i = 0, n-1 DO FOR j = 0, n-1 DO $
+;                    PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j])
+;           or equivalently, in vector notation,
+;                  PCOR = COV / (PERROR # PERROR)
 ;
 ;           If NOCOVAR is set or MPFIT terminated abnormally, then
 ;           COVAR is set to a scalar with value !VALUES.D_NAN.
@@ -275,7 +341,11 @@
 ;   DOF - number of degrees of freedom, computed as
 ;             DOF = N_ELEMENTS(DEVIATES) - NFREE
 ;         Note that this doesn't account for pegged parameters (see
-;         NPEGGED).
+;         NPEGGED).  It also does not account for data points which
+;         are assigned zero weight, for example if :
+;           * WEIGHTS[i] EQ 0, or
+;           * ERR[i] EQ infinity, or 
+;           * any of the values is "undefined" and /NAN is set.
 ;
 ;   ERRMSG - a string error or warning message is returned.
 ;
@@ -341,8 +411,13 @@
 ;                       parameter values.
 ;
 ;   MAXITER - The maximum number of iterations to perform.  If the
-;             number is exceeded, then the STATUS value is set to 5
-;             and MPFIT returns.
+;             number of calculation iterations exceeds MAXITER, then
+;             the STATUS value is set to 5 and MPFIT returns.  
+;
+;             If MAXITER EQ 0, then MPFIT does not iterate to adjust
+;             parameter values; however, the user function is evaluated
+;             and parameter errors/covariance/Jacobian are estimated
+;             before returning.
 ;             Default: 200 iterations
 ;
 ;   NAN - ignore infinite or NaN values in the Y, ERR or WEIGHTS
@@ -366,12 +441,15 @@
 ;
 ;   NPRINT - The frequency with which ITERPROC is called.  A value of
 ;            1 indicates that ITERPROC is called with every iteration,
-;            while 2 indicates every other iteration, etc.  Note that
-;            several Levenberg-Marquardt attempts can be made in a
-;            single iteration.
+;            while 2 indicates every other iteration, etc.  Be aware
+;            that several Levenberg-Marquardt attempts can be made in
+;            a single iteration.  Also, the ITERPROC is *always*
+;            called for the final iteration, regardless of the
+;            iteration number.
 ;            Default value: 1
 ;
-;   PARINFO - Provides a mechanism for more sophisticated constraints
+;   PARINFO - A one-dimensional array of structures.
+;             Provides a mechanism for more sophisticated constraints
 ;             to be placed on parameter values.  When PARINFO is not
 ;             passed, then it is assumed that all parameters are free
 ;             and unconstrained.  Values in PARINFO are never 
@@ -400,13 +478,43 @@
 ;              DOF     = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom
 ;              PCERROR = PERROR * SQRT(BESTNORM / DOF)   ; scaled uncertainties
 ;
+;   PFREE_INDEX - upon return, PFREE_INDEX contains an index array
+;                 which indicates which parameter were allowed to
+;                 vary.  I.e. of all the parameters PARMS, only
+;                 PARMS[PFREE_INDEX] were varied.
+;
+;   QUERY - if set, then MPFIT() will return immediately with one of
+;           the following values:
+;                 1 - if MIN_VERSION is not set
+;                 1 - if MIN_VERSION is set and MPFIT satisfies the minimum
+;                 0 - if MIN_VERSION is set and MPFIT does not satisfy it
+;           Default: not set.
+;
 ;   QUIET - set this keyword when no textual output should be printed
 ;           by MPFIT
 ;
-;   STATUS - an integer status code is returned.  All values other
-;            than zero can represent success.  It can have one of the
+;   STATUS - an integer status code is returned.  All values greater
+;            than zero can represent success (however STATUS EQ 5 may
+;            indicate failure to converge).  It can have one of the
 ;            following values:
 ;
+;        -18  a fatal execution error has occurred.  More information
+;             may be available in the ERRMSG string.
+;
+;        -16  a parameter or function value has become infinite or an
+;             undefined number.  This is usually a consequence of
+;             numerical overflow in the user's model function, which
+;             must be avoided.
+;
+;        -15 to -1 
+;             these are error codes that either MYFUNCT or ITERPROC
+;             may return to terminate the fitting process (see
+;             description of MPFIT_ERROR common below).  If either
+;             MYFUNCT or ITERPROC set ERROR_CODE to a negative number,
+;             then that number is returned in STATUS.  Values from -15
+;             to -1 are reserved for the user functions and will not
+;             clash with MPFIT.
+;
 ;	   0  improper input parameters.
 ;         
 ;	   1  both actual and predicted relative reductions
@@ -542,10 +650,14 @@
 ;     (the weights were not being applied), CM, 03 Sep 2007
 ;   Add COMPATIBILITY section, CM, 13 Dec 2007
 ;   Add documentation about NAN behavior, CM, 30 Mar 2009
+;   Add keywords BEST_RESIDS, CALC_FJAC, BEST_FJAC, PFREE_INDEX;
+;     update some documentation that had become stale, CM, 2010-10-28
+;   Documentation corrections, CM, 2011-08-26
+;   Additional documentation about explicit derivatives, CM, 2012-07-23
 ;
-;  $Id: mpfitfun.pro,v 1.1 2009/08/28 14:29:35 bernard Exp $
+;  $Id: mpfitfun.pro,v 1.19 2012/09/27 23:59:31 cmarkwar Exp $
 ;-
-; Copyright (C) 1997-2002, 2003, 2006, 2007, 2009, Craig Markwardt
+; Copyright (C) 1997-2002, 2003, 2006, 2007, 2009, 2010, 2011, 2012, Craig Markwardt
 ; This software is provided as is without any warranty whatsoever.
 ; Permission to use, copy, modify, and distribute modified or
 ; unmodified copies is granted, provided this copyright and disclaimer
@@ -629,6 +741,8 @@ end
  
 function mpfitfun, fcn, x, y, err, p, WEIGHTS=wts, FUNCTARGS=fa, $
                    BESTNORM=bestnorm, nfev=nfev, STATUS=status, $
+                   best_resid=best_resid, pfree_index=ifree, $
+                   calc_fjac=calc_fjac, best_fjac=best_fjac, $
                    parinfo=parinfo, query=query, CASH=cash, $
                    covar=covar, perror=perror, yfit=yfit, $
                    niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
@@ -733,11 +847,20 @@ function mpfitfun, fcn, x, y, err, p, WEIGHTS=wts, FUNCTARGS=fa, $
   result = mpfit('mpfitfun_eval', p, SCALE_FCN=scalfcn, $
                  parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$
                  covar=covar, perror=perror, $
+                 best_resid=best_resid, pfree_index=ifree, $
+                 calc_fjac=calc_fjac, best_fjac=best_fjac, $
                  niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                  ERRMSG=errmsg, quiet=quiet, _EXTRA=extra)
  
   ;; Retrieve the fit value
   yfit = temporary(mc)
+
+  ;; Rescale the Jacobian according to parameter uncertainties
+  if keyword_set(calc_fjac) AND nfree GT 0 AND status GT 0 then begin
+      ec = 1/wc  ;; Per-data-point errors (could be INF or NAN!)
+      for i = 0, nfree-1 do best_fjac[*,i] = - best_fjac[*,i] * ec
+  endif
+
   ;; Some cleanup
   xc = 0 & yc = 0 & wc = 0 & ec = 0 & mc = 0 & ac = 0
  
@@ -38,10 +38,10 @@ pro init_healpix, verbose=verbose
 healpix_sysvar = '!HEALPIX'
  
 ; Healpix version
-version = '2.11c'
+version = '2.12a'
  
 ; release data
-date = '2009-02-19'
+date = '2009-08-06'
  
 ; Healpix directory
 sv = 'HEALPIX'
-pro plotsym, psym, psize, FILL=fill,thick=thick    ;Define some plotting symbols
-;+
+pro plotsym, psym, psize, FILL=fill,thick=thick,Color = color
+   ;+
 ; NAME:
 ;     PLOTSYM
 ; PURPOSE:
@@ -8,8 +8,10 @@ pro plotsym, psym, psize, FILL=fill,thick=thick    ;Define some plotting symbols
 ;     After a symbol has been defined with PLOTSYM, a plotting command should
 ;     follow with either PSYM = 8 or !P.PSYM = 8 (see USERSYM)
 ;
+;     For additional rotationally symmetric plotting symbols, see VSYM.PRO
+;     or http://www.dfanning.com/documents/programs.html#SYMCAT
 ; CALLING SEQUENCE:
-;     PLOTSYM, PSYM,[ PSIZE, /FILL, THICK=]
+;     PLOTSYM, PSYM,[ PSIZE, /FILL, THICK=, COLOR=]
 ;
 ; INPUTS:
 ;     PSYM -  The following integer values of PSYM will create the
@@ -35,7 +37,7 @@ pro plotsym, psym, psize, FILL=fill,thick=thick    ;Define some plotting symbols
 ;             The default is 0, unfilled symbol.  Does not affect arrows
 ;             or character symbols.
 ;     THICK -  Thickness of unfilled symbols. Default is 1.
-;
+;     COLOR - Color of the symbols, Default is !P.color
 ; OUTPUTS:
 ;     None
 ;
@@ -58,12 +60,13 @@ pro plotsym, psym, psize, FILL=fill,thick=thick    ;Define some plotting symbols
 ;      Written       W. Landsman         June 1992
 ;      18-JAN-1996    Added a square symbol, HCW.
 ;      98Aug20         Added keyword thick parameter - RCB.
+;      April 2001     Added COLOR keyword    WBL
 ;-
  On_error,2
  
  if N_elements(psym) LT 1 then begin
      print,'Syntax - PLOTSYM, psym, [ size, /FILL, THICK= ]'
-     print,'  PSYM values 0 - circle, 1 - down arrow, 2 - up arrow, 3 - star
+     print,'  PSYM values 0 - circle, 1 - down arrow, 2 - up arrow, 3 - star'
      print,'       4 - triangle, 5 - upside down triangle, 6 - left arrow'
      print,'       7 - right arrow, 8 - square'
      return
@@ -90,10 +93,15 @@ pro plotsym, psym, psize, FILL=fill,thick=thick    ;Define some plotting symbols
     fill = 0
     end
 3:  begin                                     ;Star
-    r = psize
-    ang = (720. / 5*findgen(6) + 45) / !RADEG  ;Define star angles every 144 deg
-    xarr = r*cos(ang)  & yarr = r*sin(ang)
-    end
+   ang = (360. / 10 * findgen(11) + 90) / !RADEG  ;star angles every 36 deg
+   r = ang*0
+   r[2*indgen(6)] = 1.
+   cp5 = cos(!pi/5.)
+   r1 = 2. * cp5 - 1. / cp5
+   r[2*indgen(5)+1] = r1
+   r = r * psize / sqrt(!pi/4.) * 2. / (1.+r1)
+   xarr = r * cos(ang)   &   yarr = r * sin(ang)
+   end
 4:  begin                                     ;Triangle
     xarr = [-1,0,1,-1]*psize
     yarr = [-1,1,-1,-1]*psize
@@ -119,8 +127,9 @@ pro plotsym, psym, psize, FILL=fill,thick=thick    ;Define some plotting symbols
  else: message,'Unknown plotting symbol value of '+strtrim(psym,2)
  endcase
  
+ if N_elements(color) GT 0 then $
+ usersym, xarr, yarr, FILL = fill,thick=thick, color = color else $
  usersym, xarr, yarr, FILL = fill,thick=thick
-
  return
  end
  
-;
;+
; NAME:
;       STRTRANS
; PURPOSE:
;       Translate all occurences of one substring to another.
; CATEGORY:
;       text/strings
; CALLING SEQUENCE:
;       new = strtrans(oldstr,from,to,ned)
; INPUTS:
;       oldstr -- string on which to operate.              in
;                 May be an array.
;       from   -- substrings to be translated. May be      in
;                 an array.
;       to     -- what strings in from should be           in
;                 translated to. May be an array.
; KEYWORD PARAMETERS:
;       /HELP  -- Set this to print useful message and 
;                 exit.
; OUTPUTS:
;       new    -- Translated string. Array if oldstr is    out          
;                 an array.
;       ned    -- number of substitutions performed in     out
;                 oldstr.  Array if oldstr is an array.
; COMMON BLOCKS:
; SIDE EFFECTS:
; NOTES:
;       - Any of old, from, and to can be arrays.  
;       - from and to must have the same number of elements.
; EXAMPLE:
;       inp='Many*bad!chars+in_here'
;       from=['*','!','+','_']
;       to  =[' ',' ',' ',' ']
;       out = strtrans(inp,from,to,ned)
;       Will produce out='Many bad chars in here', and set ned to 4.
; MODIFICATION HISTORY:
;       $Id: strtrans.pro,v 1.4 2001/11/21 19:13:23 mcraig Exp $
;       $Log: strtrans.pro,v $
;       Revision 1.4  2001/11/21 19:13:23  mcraig
;       Changed str_sep to strsplit. The former is now considered obsolete by RSI.
;
;       Revision 1.3  1996/06/14 20:00:27  mcraig
;       Updated Copyright info.
;
;       Revision 1.2  1996/05/09 00:22:17  mcraig
;       Sped up significantly by using str_sep to handle the translation.  No longer
;       relies on routines fromother user libraries.
;
;       Revision 1.1  1996/01/31 18:47:37  mcraig
;       Initial revision
;
; RELEASE:
;       $Name: Rel_2_1 $
;
; COPYRIGHT:
;  Copyright (C) 1996 The Regents of the University of California, All
;  Rights Reserved.  Written by Matthew W. Craig.
;  See the file COPYRIGHT for restrictions on distrubting this code.
;  This code comes with absolutely NO warranty; see DISCLAIMER for details.
;-
;
FUNCTION strtrans, InputString, from, to, ned,  $
                   HELP=Help

; Bomb out to caller if error.
    On_error, 2

; Offer help if we don't have at least InputString, from, and to, or
; if the user asks for it.
    IF (n_params() LT 3) OR keyword_set(help) THEN BEGIN
        offset = '   '
        print, offset+'Translate all occurences of one substring to another.'
        print, offset+'new = strtrans(oldstr,from,to,ned)'
        print, offset+'Inputs:'
        print, offset+offset+'oldstr -- string on which to operate.              in'
        print, offset+offset+'          May be an array.'
        print, offset+offset+'from   -- substrings to be translated. May be      in'
        print, offset+offset+'          an array.'
        print, offset+offset+'to     -- what strings in from should be           in'
        print, offset+offset+'          translated to. May be an array.'
        print, offset+'Outputs:'
        print, offset+offset+'new    -- Translated string. Array if oldstr is    out'
        print, offset+offset+'          an array.'
        print, offset+offset+'ned    -- number of substitutions performed in     out'
        print, offset+offset+'          oldstr.  Array if oldstr is an array.'
        print, offset+'Notes:'
        print, offset+offset+'- Any of old, from, and to can be arrays. ' 
        print, offset+offset+'- from and to must have the same number of elements.'
        return, -1
    ENDIF 
    
    strn = InputString

;  Check that From/To have same number of elements.  RETURN if they don't.
    NFrom = n_elements(from)
    NTo = n_elements(to)
    IF (NFrom EQ 0) OR (NTo EQ 0) THEN return, strn
    IF NFrom NE NTo THEN BEGIN
        print,'Error: Number of elements in from/to unequal'
        return,-1
    ENDIF

;  Make sure there are no null strings in From.  RETURN if there are.   
    FromLen = strlen(From)
    IF (total(FromLen EQ 0) GT 0) THEN BEGIN
        print, 'Error: elements of From must have nonzero length.'
        return, -1
    ENDIF 

    NStrings = n_elements(strn)
    ned = lonarr(NStrings)
    tmpned = 0L
; Say strn='a#b#c', from='#' and to='@'.  Then the approach here is to
; first split strn at all occurances of '#', then recombine the pieces
; with '@' inserted instead.  Do this for all elements of strn, and
; all elements of from.
    FOR i = 0L, NStrings-1 DO BEGIN
        ned(i) = 0L
        st=strn(i)
        FOR j=0L, NFrom-1 DO BEGIN
            strn[i]=textoidl_str_replace(strn[i],from(j),to(j))       ;jpb patch using JD's routine
;            SepStr = strsplit(strn(i), from(j))             ;jpb patch
;            pos=strposmulti(strn(i),from(j),count)
;            SepStr = strsplit(strn(i), from(j),/extract)             ;jpb patch
;            NSubs = n_elements(SepStr) - 1
;            IF count NE 0 THEN BEGIN
;              stop
;              SepStr = strsplit(strn(i), from(j),/extract)
;              NSubs = n_elements(SepStr) - 1
;              strn(i) = SepStr(0)
;              FOR k=1L, NSubs DO strn(i) = strn(i) + To(j) + SepStr(k)
;              ned(i) =  ned(i) + NSubs
;            ENDIF
        ENDFOR 
    ENDFOR

    return, strn
END 
 \ No newline at end of file
+;
;+
; NAME:
;       STRTRANS
; PURPOSE:
;       Translate all occurences of one substring to another.
; CATEGORY:
;       text/strings
; CALLING SEQUENCE:
;       new = strtrans(oldstr,from,to,ned)
; INPUTS:
;       oldstr -- string on which to operate.              in
;                 May be an array.
;       from   -- substrings to be translated. May be      in
;                 an array.
;       to     -- what strings in from should be           in
;                 translated to. May be an array.
; KEYWORD PARAMETERS:
;       /HELP  -- Set this to print useful message and 
;                 exit.
; OUTPUTS:
;       new    -- Translated string. Array if oldstr is    out          
;                 an array.
;       ned    -- number of substitutions performed in     out
;                 oldstr.  Array if oldstr is an array.
; COMMON BLOCKS:
; SIDE EFFECTS:
; NOTES:
;       - Any of old, from, and to can be arrays.  
;       - from and to must have the same number of elements.
; EXAMPLE:
;       inp='Many*bad!chars+in_here'
;       from=['*','!','+','_']
;       to  =[' ',' ',' ',' ']
;       out = strtrans(inp,from,to,ned)
;       Will produce out='Many bad chars in here', and set ned to 4.
; MODIFICATION HISTORY:
;       $Id: strtrans.pro,v 1.4 2001/11/21 19:13:23 mcraig Exp $
;       $Log: strtrans.pro,v $
;       Revision 1.4  2001/11/21 19:13:23  mcraig
;       Changed str_sep to strsplit. The former is now considered obsolete by RSI.
;
;       Revision 1.3  1996/06/14 20:00:27  mcraig
;       Updated Copyright info.
;
;       Revision 1.2  1996/05/09 00:22:17  mcraig
;       Sped up significantly by using str_sep to handle the translation.  No longer
;       relies on routines fromother user libraries.
;
;       Revision 1.1  1996/01/31 18:47:37  mcraig
;       Initial revision
;
; RELEASE:
;       $Name: Rel_2_1 $
;
; COPYRIGHT:
;  Copyright (C) 1996 The Regents of the University of California, All
;  Rights Reserved.  Written by Matthew W. Craig.
;  See the file COPYRIGHT for restrictions on distrubting this code.
;  This code comes with absolutely NO warranty; see DISCLAIMER for details.
;-
;
FUNCTION strtrans, InputString, from, to, ned,  $
                   HELP=Help

; Bomb out to caller if error.
    On_error, 2

; Offer help if we don't have at least InputString, from, and to, or
; if the user asks for it.
    IF (n_params() LT 3) OR keyword_set(help) THEN BEGIN
        offset = '   '
        print, offset+'Translate all occurences of one substring to another.'
        print, offset+'new = strtrans(oldstr,from,to,ned)'
        print, offset+'Inputs:'
        print, offset+offset+'oldstr -- string on which to operate.              in'
        print, offset+offset+'          May be an array.'
        print, offset+offset+'from   -- substrings to be translated. May be      in'
        print, offset+offset+'          an array.'
        print, offset+offset+'to     -- what strings in from should be           in'
        print, offset+offset+'          translated to. May be an array.'
        print, offset+'Outputs:'
        print, offset+offset+'new    -- Translated string. Array if oldstr is    out'
        print, offset+offset+'          an array.'
        print, offset+offset+'ned    -- number of substitutions performed in     out'
        print, offset+offset+'          oldstr.  Array if oldstr is an array.'
        print, offset+'Notes:'
        print, offset+offset+'- Any of old, from, and to can be arrays. ' 
        print, offset+offset+'- from and to must have the same number of elements.'
        return, -1
    ENDIF 
    
    strn = InputString

;  Check that From/To have same number of elements.  RETURN if they don't.
    NFrom = n_elements(from)
    NTo = n_elements(to)
    IF (NFrom EQ 0) OR (NTo EQ 0) THEN return, strn
    IF NFrom NE NTo THEN BEGIN
        print,'Error: Number of elements in from/to unequal'
        return,-1
    ENDIF

;  Make sure there are no null strings in From.  RETURN if there are.   
    FromLen = strlen(From)
    IF (total(FromLen EQ 0) GT 0) THEN BEGIN
        print, 'Error: elements of From must have nonzero length.'
        return, -1
    ENDIF 

    NStrings = n_elements(strn)
    ned = lonarr(NStrings)
    tmpned = 0L
; Say strn='a#b#c', from='#' and to='@'.  Then the approach here is to
; first split strn at all occurances of '#', then recombine the pieces
; with '@' inserted instead.  Do this for all elements of strn, and
; all elements of from.
    FOR i = 0L, NStrings-1 DO BEGIN
        ned(i) = 0L
        st=strn(i)
        FOR j=0L, NFrom-1 DO BEGIN
            strn(i)=textoidl_str_replace(strn(i),from(j),to(j))       ;jpb patch using JD's routine
;            SepStr = strsplit(strn(i), from(j))             ;jpb patch
;            pos=strposmulti(strn(i),from(j),count)
;            SepStr = strsplit(strn(i), from(j),/extract)             ;jpb patch
;            NSubs = n_elements(SepStr) - 1
;            IF count NE 0 THEN BEGIN
;              stop
;              SepStr = strsplit(strn(i), from(j),/extract)
;              NSubs = n_elements(SepStr) - 1
;              strn(i) = SepStr(0)
;              FOR k=1L, NSubs DO strn(i) = strn(i) + To(j) + SepStr(k)
;              ned(i) =  ned(i) + NSubs
;            ENDIF
        ENDFOR 
    ENDFOR

    return, strn
END 
 \ No newline at end of file
-;
;+
; NAME:
;       TEXTOIDL
; PURPOSE:
;       Convert a valid TeX string to a valid IDL string for plot labels.
; CATEGORY:
;       text/strings
; CALLING SEQUENCE:
;       new = textoidl(old)
; INPUTS:
;       old            -- TeX string to be converted.  Will not be     in
;                         modified.  old may be a string array.
; KEYWORD PARAMETERS:
;       FONT           -- Set to 0 to use hardware font, -1 to use 
;                         vector.  Note that the only hardware font 
;                         supported is PostScript.
;       /TEX_SEQUENCES -- return the available TeX sequences
;       /HELP          -- print out info on use of the function
;                         and exit.
; OUTPUTS:
;       new            -- IDL string corresponding to old.             out
; COMMON BLOCKS:
; SIDE EFFECTS:
; NOTES:
;       - Use the procedure SHOWTEX to get a list of the available TeX
;         control sequences.  
;       - The only hardware font for which translation is available is
;         PostScript. 
;       - The only device for which hardware font'
;         translation is available is PostScript.'
;       - The FONT keyword overrides the font selected'
;         by !p.font'
; EXAMPLE:
;       out = TeXtoIDL('\Gamma^2 + 5N_{ed}')
;       The string out may be used in XYOUTS or other IDL text
;       display routines.  It will be an uppercase Gamma, with an
;       exponent of 2, then a plus sign, then an N with the subscript
;       ed.
; MODIFICATION HISTORY:
;       $Id: textoidl.pro,v 1.4 1996/06/14 20:00:27 mcraig Exp $
;       $Log: textoidl.pro,v $
;       Revision 1.4  1996/06/14 20:00:27  mcraig
;       Updated Copyright info.
;
;       Revision 1.3  1996/05/09 00:22:17  mcraig
;       Added error handling, cleaned up documentation.
;
;       Revision 1.2  1996/02/08 18:52:50  mcraig
;       Added ability to use hardware fonts for PostScript device.
;
;       Revision 1.1  1996/01/31 18:47:37  mcraig
;       Initial revision
;
; RELEASE:
;       $Name: Rel_2_1 $
;
; COPYRIGHT:
;  Copyright (C) 1996 The Regents of the University of California, All
;  Rights Reserved.  Written by Matthew W. Craig.
;  See the file COPYRIGHT for restrictions on distrubting this code.
;  This code comes with absolutely NO warranty; see DISCLAIMER for details.
;-
;
FUNCTION Textoidl, InputString, $
                   FONT=fnt, $
                   HELP=hlp, $
                   TEX_SEQUENCES=tex_seq

;  Return to caller if there is an error.
    On_error, 2
;  We begin by deciding on the font.  PostScript = 0 means use vector.
    PostScript = 0
    IF n_elements(fnt) EQ 0 THEN BEGIN     ; get font from !p.font
        IF !p.font NE -1 THEN BEGIN        ; User wants hardware font.
            PostScript=1
        ENDIF
    ENDIF ELSE BEGIN                       ; get font from FONT keyword
        IF fnt NE -1 THEN PostScript = 1
    ENDELSE

;  Bomb out if user wants non-PostScript hardware font.
    IF (PostScript EQ 1) AND (!d.name NE 'PS') THEN BEGIN   
                                              ; Device isn't postscript 
                                              ; and user wants hardware
                                              ; font.  Not good.
        print,'Warning: No translation for device: ',!d.name
        return,InputString               
    ENDIF 
    
    IF keyword_set (tex_seq) THEN BEGIN
        table=textable()
        return,table(0,*)
    ENDIF 

    IF keyword_set(hlp) OR (n_params() EQ 0) THEN BEGIN
        print, '   Convert a TeX string to an IDL string'
        print, '   new = TeXtoIDL(old)'
        print, '     old = TeX string to translate.                 in'
        print, '     new = resulting IDL string.                    out'
        print, '   Keywords:'
        print, '      FONT       set to -1 to translate for vector fonts '
        print, '                 (DEFAULT) .  Set to 0 to translate for'
        print, '                 hardware font.'
        print, '      /TEX_SEQUENCES -- return the available TeX sequences'
        print, '      /HELP      print this message and exit.'
        print, '   NOTES:  '
        print, '      - Use SHOWTEX to obtain a list of the available'
        print, '        TeX control sequences.'
        print, '      - old may be a string array.  If so, new is too.'
        print, '      - The only device for which hardware font'
        print, '        translation is available is PostScript.'
        print, '      - The FONT keyword overrides the font selected'
        print, '        by !p.font'
        return, -1
    ENDIF
    
; PostScript has been set to 1 if PostScript fonts are desired.
    strn = InputString
    table = textable(POSTSCRIPT=PostScript)
;   Greek sub/superscripts need to be protected by putting braces
;   around them if they are unbraced.  This will have the result the
;   it will be difficult to use \ as a sub/superscript.  Get over it.
    strn =  strtrans(strn, '^'+table(0, *), '^{'+table(0, *)+'}')
    strn =  strtrans(strn, '_'+table(0, *), '_{'+table(0, *)+'}')

;  First we translate Greek letters and the like.  This makes guessing
;  alignment of sub/superscripts easier, as all special characters will then
;  be one character long.
    strn = strtrans(strn, table(0, *), table(1, *))

    FOR i = 0L, n_elements(strn)-1 DO $
      strn[i] = translate_sub_super(strn[i]) ; Take care of sub/superscripts

    return,strn
END 
 \ No newline at end of file
+;
;+
; NAME:
;       TEXTOIDL
; PURPOSE:
;       Convert a valid TeX string to a valid IDL string for plot labels.
; CATEGORY:
;       text/strings
; CALLING SEQUENCE:
;       new = textoidl(old)
; INPUTS:
;       old            -- TeX string to be converted.  Will not be     in
;                         modified.  old may be a string array.
; KEYWORD PARAMETERS:
;       FONT           -- Set to 0 to use hardware font, -1 to use 
;                         vector.  Note that the only hardware font 
;                         supported is PostScript.
;       /TEX_SEQUENCES -- return the available TeX sequences
;       /HELP          -- print out info on use of the function
;                         and exit.
; OUTPUTS:
;       new            -- IDL string corresponding to old.             out
; COMMON BLOCKS:
; SIDE EFFECTS:
; NOTES:
;       - Use the procedure SHOWTEX to get a list of the available TeX
;         control sequences.  
;       - The only hardware font for which translation is available is
;         PostScript. 
;       - The only device for which hardware font'
;         translation is available is PostScript.'
;       - The FONT keyword overrides the font selected'
;         by !p.font'
; EXAMPLE:
;       out = TeXtoIDL('\Gamma^2 + 5N_{ed}')
;       The string out may be used in XYOUTS or other IDL text
;       display routines.  It will be an uppercase Gamma, with an
;       exponent of 2, then a plus sign, then an N with the subscript
;       ed.
; MODIFICATION HISTORY:
;       $Id: textoidl.pro,v 1.4 1996/06/14 20:00:27 mcraig Exp $
;       $Log: textoidl.pro,v $
;       Revision 1.4  1996/06/14 20:00:27  mcraig
;       Updated Copyright info.
;
;       Revision 1.3  1996/05/09 00:22:17  mcraig
;       Added error handling, cleaned up documentation.
;
;       Revision 1.2  1996/02/08 18:52:50  mcraig
;       Added ability to use hardware fonts for PostScript device.
;
;       Revision 1.1  1996/01/31 18:47:37  mcraig
;       Initial revision
;
; RELEASE:
;       $Name: Rel_2_1 $
;
; COPYRIGHT:
;  Copyright (C) 1996 The Regents of the University of California, All
;  Rights Reserved.  Written by Matthew W. Craig.
;  See the file COPYRIGHT for restrictions on distrubting this code.
;  This code comes with absolutely NO warranty; see DISCLAIMER for details.
;-
;
FUNCTION Textoidl, InputString, $
                   FONT=fnt, $
                   HELP=hlp, $
                   TEX_SEQUENCES=tex_seq

;  Return to caller if there is an error.
    On_error, 2
;  We begin by deciding on the font.  PostScript = 0 means use vector.
    PostScript = 0
    IF n_elements(fnt) EQ 0 THEN BEGIN     ; get font from !p.font
        IF !p.font NE -1 THEN BEGIN        ; User wants hardware font.
            PostScript=1
        ENDIF
    ENDIF ELSE BEGIN                       ; get font from FONT keyword
        IF fnt NE -1 THEN PostScript = 1
    ENDELSE

;  Bomb out if user wants non-PostScript hardware font.
    IF (PostScript EQ 1) AND (!d.name NE 'PS') THEN BEGIN   
                                              ; Device isn't postscript 
                                              ; and user wants hardware
                                              ; font.  Not good.
        print,'Warning: No translation for device: ',!d.name
        return,InputString               
    ENDIF 
    
    IF keyword_set (tex_seq) THEN BEGIN
        table=textable()
        return,table(0,*)
    ENDIF 

    IF keyword_set(hlp) OR (n_params() EQ 0) THEN BEGIN
        print, '   Convert a TeX string to an IDL string'
        print, '   new = TeXtoIDL(old)'
        print, '     old = TeX string to translate.                 in'
        print, '     new = resulting IDL string.                    out'
        print, '   Keywords:'
        print, '      FONT       set to -1 to translate for vector fonts '
        print, '                 (DEFAULT) .  Set to 0 to translate for'
        print, '                 hardware font.'
        print, '      /TEX_SEQUENCES -- return the available TeX sequences'
        print, '      /HELP      print this message and exit.'
        print, '   NOTES:  '
        print, '      - Use SHOWTEX to obtain a list of the available'
        print, '        TeX control sequences.'
        print, '      - old may be a string array.  If so, new is too.'
        print, '      - The only device for which hardware font'
        print, '        translation is available is PostScript.'
        print, '      - The FONT keyword overrides the font selected'
        print, '        by !p.font'
        return, -1
    ENDIF
    
; PostScript has been set to 1 if PostScript fonts are desired.
    strn = InputString
    table = textable(POSTSCRIPT=PostScript)
;   Greek sub/superscripts need to be protected by putting braces
;   around them if they are unbraced.  This will have the result the
;   it will be difficult to use \ as a sub/superscript.  Get over it.
    strn =  strtrans(strn, '^'+table(0, *), '^{'+table(0, *)+'}')
    strn =  strtrans(strn, '_'+table(0, *), '_{'+table(0, *)+'}')

;  First we translate Greek letters and the like.  This makes guessing
;  alignment of sub/superscripts easier, as all special characters will then
;  be one character long.
    strn = strtrans(strn, table(0, *), table(1, *))

    FOR i = 0L, n_elements(strn)-1 DO $
      strn(i) = translate_sub_super(strn(i)) ; Take care of sub/superscripts

    return,strn
END 
 \ No newline at end of file