Merge branch 'master' of https://gitlab.irap.omp.eu/OV-GSO-DC/dustem-wrapper_idl

Annie Hughes
2 parents 9c053027 933ab1ab
Showing 8 changed files with 54 additions and 7 deletions Show diff stats
instrument_description.xcat
src/idl/Markwardt_for_Dustemwrap/dustem_mpfit.pro
src/idl/dustem_init_parinfo.pro
src/idl/dustem_mpfit_data.pro
src/idl/dustem_mpfit_run.pro
src/idl/dustem_plugin_stellar_population.pro
src/idl/dustem_read_filters.pro
src/idl/dustem_show_file_dependencies.pro
@@ -174,11 +174,11 @@ SPASS      SPASS1     130174.74  MJy/sr      NULL         0.14833333          5.
 \========= ZOA-FF (free free from Alves, M. I. R., et al, 2015, MNRAS, 450, 2025) ==========
 \information taken from Alves, M. I. R., et al, 2015, MNRAS, 450, 2025
 \flux convention does not really matter (narrow band of ? MHz for 1.4 GHz)
-CHIPASS    ZOA_FF1    214286.    MJy/sr      NULL         0.24                5.0
+\CHIPASS    ZOA_FF1    214286.    MJy/sr      NULL         0.24                5.0
 \========= CHI-PASS ==========
 \information taken from Calabretta, M. et al. 2014, PASA, 31, 7
 \flux convention does not really matter (narrow band of ? MHz for 1.4 GHz)
-CHIPASS    CHIPASS1   214286.     MJy/sr      NULL        0.24                5.0
+CHIPASS    CHIPASS1   214286.    MJy/sr      NULL        0.24                5.0
 \========= IRAC ====================================================================================================
 \resolution taken to be 1.7, 1.7, 1.9, 2" (eg Meixner et al)
 IRAC       IRAC1      3.55000    MJy/sr      IRAC         0.000472222
@@ -3571,7 +3571,19 @@ function dustem_mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
       ;#iterations
       !dustem_iter.act=iter
       ;#current parameter values
-      (*!dustem_fit).current_param_values=ptr_new(x*(*(*!dustem_fit).param_init_values));So that parameters aren't just passed directly from MPFIT. ALSO allows to test with !dustem_current
+      ;JPB: Below is because whene there are fixed parameters, x and (*!dustem_fit) don't have same number of elements.
+      ;JPB: in principle, we should not have modified the code of mpfit for dustem.
+      ;JPB: The same should be achievable using keywords as iterproc
+      IF keyword_set(parinfo) THEN BEGIN
+      	xx=dblarr(n_elements(parinfo))
+      	ind1=where(parinfo.fixed EQ 0,count1,complement=ind2,Ncomplement=count2)
+      	IF count1 NE 0 THEN xx[ind1]=x
+      	IF count2 NE 0 THEN xx[ind2]=parinfo[ind2].value
+      ENDIF ELSE BEGIN
+      	xx=x
+      ENDELSE
+      ;stop
+      (*!dustem_fit).current_param_values=ptr_new(xx*(*(*!dustem_fit).param_init_values));So that parameters aren't just passed directly from MPFIT. ALSO allows to test with !dustem_current
 ;     ;#Current parameter errors
  
   if n_elements(fnorm) GT 0 AND n_elements(fnorm1) GT 0 then begin
@@ -7,7 +7,7 @@ PRO dustem_init_parinfo,pd,iv, $
  
 (*!dustem_fit).param_init_values=ptr_new(iv)
  
-(*!dustem_fit).current_param_values = ptr_new(iv/iv)    
+(*!dustem_fit).current_param_values = ptr_new(iv/iv)  ;These are normalized to the initial values  
 ;stop
 dustem_set_func_ind,pd,iv
  
@@ -132,6 +132,8 @@ p_min = dustem_mpfitfun(func_name,wav,values,sigma, $
                  maxiter=use_Nitermax,gtol=gtol,errmsg=errmsg,functargs=_extra,nocatch=nocatch)
  
 ;stop                 
+;!dustem_iter.act=niter_completed
+
 message,errmsg,/info
 message,'mpfitfun Status='+strtrim(status,2),/info
 p_dim=p_min*(*(*!dustem_fit).param_init_values)
@@ -151,7 +153,7 @@ IF !dustem_noobj THEN BEGIN
 	dustemwrap_plot,p_dim,stp,_extra=_extra
 ENDELSE
  
-     the_end:
+the_end:
 RETURN, p_dim
  
 END
@@ -111,8 +111,10 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
  
 		'SED'	: BEGIN 
 	        dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec)
+          ;help,*(*!dustem_fit).CURRENT_PARAM_VALUES
+          ;stop
 	        chi2_sed   = total(((*(*!dustem_data).sed).values-dustem_sed)^2/(*(*!dustem_data).sed).sigma^2)
-			n_sed = n_elements((*(*!dustem_data).sed).values)
+			    n_sed = n_elements((*(*!dustem_data).sed).values)
 	       	rchi2_sed  = chi2_sed / n_sed
 	       	wrchi2_sed  = rchi2_sed * (*!fit_rchi2_weight).sed
 	       	message,"chi2 SED = "+strtrim(chi2_sed,2)+"  rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed
@@ -92,10 +92,14 @@ IF keyword_set(key) THEN BEGIN
  
     ;JPB: This is highly non standard and makes Fawlty go "Segmentation fault"
     ;lumclass = (strmid(key,2)).remove(-1) ; extracting the luminosity class strings 
+    ;IC: Yes but we need to take into account when luminosity class is two characters and not just one. ie: IV instead of V. 
+    ;But since we only have MS (V) data, (lumclass[i]=strmid(key[i],2,1)) should be ok but I'm correcting it for future input.
     lumclass=strarr(n_elements(key))
     FOR i=0L,n_elements(key)-1 DO BEGIN
-        lumclass[i]=strmid(key[i],2,1)
+        lumclass[i]=strmid(key[i],2,strlen(key[i])-3)
     ENDFOR
+    
+    
     params = strmid(key,0,/reverse_offset) ; extracting the indices (in the 'key' array) of the stellar parameters to fit  
  
     ; Locating the O3 stellar population parameters - this new code structure accounts for the presence of different luminosity classes.
@@ -816,6 +816,30 @@ spass_struct={Name:&#39;SPASS&#39;,Nbands:Nband,filter_names:filt_names,central_waveleng
             use_frequencies:use_nu,use_transmissions:use_T,use_wavelengths:use_w,use_wmin:dblarr(Nband),use_wmax:dblarr(Nband),cc_method:'dustem_cc_spass'}
  
 ;======================================================
+;==== CHIPASS
+;======================================================
+IF print_messages THEN message,'Reading CHIPASS filters',/continue
+
+Nband=1
+filt_names=['CHIPASS1']
+wfilt_chipass=dustem_filter2wav(filt_names)
+wav_ptr=ptrarr(Nband)
+nu_ptr=ptrarr(Nband) & T_ptr=ptrarr(Nband) & wav_ptr=ptrarr(Nband)
+use_nu=ptrarr(Nband) & use_T=ptrarr(Nband) & use_w=ptrarr(Nband)
+
+;===== This is set on purpose to the SPASS directory.
+st=read_xcat(dir_spass+'SPASS_transmission.xcat',/silent)
+
+freq=st.frequency*1.e9  ;Hz
+Trans=st.transmission
+
+T_ptr(0)=ptr_new(Trans/max(Trans)) & nu_ptr(0)=ptr_new(freq) & wav_ptr(0)=ptr_new(cmic/freq)  
+
+chipass_struct={Name:'CHIPASS',Nbands:Nband,filter_names:filt_names,central_wavelengths:wfilt_chipass,central_frequencies:cmic/wfilt_spass, $
+            filter_wavelengths:wav_ptr,filter_frequencies:nu_ptr,filter_transmissions:T_ptr, $
+            use_frequencies:use_nu,use_transmissions:use_T,use_wavelengths:use_w,use_wmin:dblarr(Nband),use_wmax:dblarr(Nband),cc_method:'dustem_cc_spass'}
+
+;======================================================
 ;==== NIKA2
 ;======================================================
 IF print_messages THEN message,'Reading NIKA2 filters',/continue
@@ -1045,6 +1069,7 @@ dm_filter_struct={ $
                  LABOCA:laboca_struct, $
                  GISMO:gismo_struct, $
                  SPASS:spass_struct, $
+                 CHIPASS:chipass_struct, $
                  NIKA2:nika2_struct, $
                  SCUBA2:scuba2_struct, $
                  IRS:irs_struct, $
@@ -72,7 +72,9 @@ dustem_fit_ext_pol_example,Nitermax=1,fits_save_and_restore=&#39;/tmp/mysavefile3.fi
 dustem_fit_sed_ext_pol_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile4.fits'
 dustem_fit_intensity_mbb_example,Nitermax=1,postscript='/tmp/dustemwrap_postcript_example.ps'
 ;=== for some reason plot goes red here under fawlty
+;=== and fails when run with only the minimal dustemwrap installed (stellar . txt file is missing)
 dustem_myisrf_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile5.fits'
+;=== fails when run with only the minimal dustemwrap installed (stellar . txt file is missing)
 dustem_stellarpopisrf_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile6.fits'
 dustem_fit_intensity_example,Nitermax=1,/noobj