OV-GSO-DC / dustem-wrapper_idl

Commit c76766dd80c32a1b35799526a7f247b35753655f

Authored by rmaris
1 parent 2a89c8d9
Exists in master

First commit

Showing 1 changed file with 143 additions and 0 deletions   Show diff stats
Inline Side-by-side
LabTools/IRAP/RM/phangs_fits_seds_sizedist.pro 0 → 100644
  Diff comments   View file @c76766d
... ... @@ -0,0 +1,143 @@
  1 +PRO phangs_fits_seds_sizedist,show_sed=show_sed $
  2 + ,normalize=normalize $
  3 + ,weighted=weighted
  4 +
  5 +;RUN EXAMPLE : phangs_fits_seds_sizedist, /show_sed, /norma,/weighted
  6 +
  7 +dustem_define_la_common
  8 +
  9 +dir=!phangs_data_dir+'/ISRF/GRIDS/'
  10 +dir_to_fit = '/home/rmaris/Desktop/Soft_Libraries/dustem-wrapper_idl/LabTools/IRAP/RM/tmp/sed_size.fits'
  11 +
  12 +isrf_class = 0
  13 +isrf_class_str='_isrfclass'+strtrim(isrf_class,2)
  14 +model='DL07'
  15 +;table_name=dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+isrf_class_str+'.fits'
  16 +table_name=dir+'DL07_logn_MuseISRF_JWST_adaptedG0_Ypah1added_Yvsgadded_4Phangs_noionis'+isrf_class_str+'.fits'
  17 +table_to_fit = dir_to_fit
  18 +; table_name = dir + 'DBP90_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs_noionis'+isrf_class_str+'.fits'
  19 +stp=mrdfits(table_to_fit ,1,hp)
  20 +st=mrdfits(table_to_fit ,2,h)
  21 +Nst=n_elements(st)
  22 +Nseds = n_elements(st.(0))
  23 +Nparams=4
  24 +Ntrueparams=2
  25 +predicted_params=fltarr(Nseds,Nparams)
  26 +true_params=fltarr(Nseds,Ntrueparams)
  27 +rchi2_res = fltarr(Nseds)
  28 +
  29 +use_filters_names=['F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W','PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW']
  30 +filters=dustem_filter_names2filters(use_filters_names)
  31 +Nfilters=n_elements(filters)
  32 +
  33 +array_seds = fltarr(Nseds,Nfilters)
  34 +
  35 +FOR j=0L,Nseds-1 DO BEGIN
  36 +
  37 + sed=dustem_initialize_sed(Nfilters)
  38 + ;sed=fltarr(Nfilters)
  39 + tags=tag_names(st)
  40 + FOR i=0L,Nfilters-1 DO BEGIN
  41 + ind=where(tags EQ 'I'+filters[i],count)
  42 + sed[i].filter=filters[i]
  43 + sed[i].stokesI=st[j].(ind[0])
  44 + ENDFOR
  45 + sed.sigmaii=(sed.stokesI*0.1)^2
  46 + FOR i=0L,Ntrueparams-1 DO BEGIN
  47 + true_params[j,i]=stp[j].(i)
  48 + ENDFOR
  49 + vsg_id=2L
  50 + vsg_id_str=strtrim(vsg_id,2)
  51 + fixed_param_desc= ['dustem_plugin_phangs_class_isrf_2'];,'(*!dustem_params).grains['+vsg_id_str+'].MDUST_O_MH']
  52 + fixed_param_val= [1.];,0.0016600000]
  53 +
  54 + array_seds[j,*]=sed.stokesI
  55 +
  56 + params=dustem_brute_force_fit(sed,table_name,filters,fact=fact,chi2=chi2,normalize=normalize,rchi2=rchi2,show_sed=show_sed $
  57 + ,params_hit=params_hit,params_uncertainties=params_uncertainties,params_min=params_min,params_max=params_max $
  58 + ,fixed_parameters_description=fixed_param_desc,fixed_parameters_values=fixed_param_val $
  59 + ,weighted_params=weighted_params,reset=reset,best_sed=best_sed,best_grid_sed=best_grid_sed,status=status,weightning=weightning $
  60 + ,weighted_sed=weighted_sed,weighted_chi2=weighted_chi2,rweighted_chi2=rweighted_chi2,marginalize_method=marginalize_method)
  61 +
  62 + ;print,params,fact
  63 + ;IF fact GT 2 THEN BEGIN
  64 + ; print,ii,fact,chi2
  65 + ; stop
  66 + ;ENDIF
  67 +
  68 + print,' true_params = ', true_params[j,*]
  69 + print,'fitted_params= ', params
  70 + print,'rweighted_chi2= ', rchi2
  71 + ;print,'best_fit_seds = ', best_sed
  72 + print,'============================================================================='
  73 +
  74 +
  75 +
  76 + FOR i=0L,Nparams-1 DO BEGIN
  77 + predicted_params[j,*]=weighted_params
  78 + ENDFOR
  79 +
  80 + rchi2_res[j] = rchi2
  81 +
  82 +
  83 +ENDFOR
  84 +
  85 +Ypah0 = predicted_params[*,1]
  86 +Ypah1 = predicted_params[*,2]
  87 +Ypah = Ypah0 + Ypah1
  88 +Xion =Ypah1/Ypah
  89 +amax_Ypah0 = true_params[*,0]
  90 +amax_Ypah1 = true_params[*,1]
  91 +
  92 +
  93 +
  94 +amax = amax_Ypah0 + amax_Ypah1
  95 +
  96 +ind = where(amax_Ypah0 EQ amax_Ypah1,count)
  97 +
  98 +
  99 +xsize=500
  100 +ysize=500
  101 +!P.Position = [1., 1., 1., 1.]
  102 +
  103 +window, 0 ,xsize=500,ysize=500
  104 +cgplot,amax,Xion,psym='Filled Circle',color='blue',/ysty,xtit='amax',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)]
  105 +cgplot,amax[ind],Xion[ind],psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)],/overplot
  106 +
  107 +ind=where(amax_Ypah1 EQ amax_Ypah1[0] ,count)
  108 +window, 1,xsize=500,ysize=500
  109 +cgplot,amax_Ypah0[ind],Xion,psym='Filled Circle',color='blue',/ysty,xtit='amax_Ypah0',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)]
  110 +
  111 +ind=where(amax_Ypah0 EQ amax_Ypah0[0] ,count)
  112 +window, 2,xsize=500,ysize=500
  113 +cgplot,amax_Ypah1[ind],Xion,psym='Filled Circle',color='blue',/ysty,xtit='amax_Ypah1',ytit='Xion',/xlog,yrange=[la_min(Xion),la_max(Xion)]
  114 +
  115 +window, 3,xsize=500,ysize=500
  116 +cgplot,amax,Ypah1,psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='Ypah1',/xlog ,yrange=[la_min(Ypah1),la_max(Ypah1)]
  117 +
  118 +window, 4,xsize=500,ysize=500
  119 +cgplot,amax,Ypah0,psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='Ypah0',/xlog ,yrange=[la_min(Ypah0),la_max(Ypah0)]
  120 +
  121 +window, 5,xsize=500,ysize=500
  122 +cgplot,amax,rchi2_res,psym='Filled Circle',color='black',xtit='amax',ytit='rchi2',/xlog,yrange=[la_min(rchi2_res),la_max(rchi2_res)]
  123 +cgplot,amax[ind],rchi2_res[ind],psym='Filled Circle',color='red',/ysty,xtit='amax',ytit='rchi2',/xlog,yrange=[la_min(rchi2_res),la_max(rchi2_res)],/overplot
  124 +
  125 +window, 6,xsize=500,ysize=500
  126 +cgplot,amax,array_seds[*,1],psym='Filled Circle',color='black',xtit='amax',ytit='F335M',/xlog,yrange=[la_min(array_seds[*,1]),la_max(array_seds[*,1])]
  127 +
  128 +window, 7,xsize=500,ysize=500
  129 +cgplot,amax,array_seds[*,3],psym='Filled Circle',color='black',xtit='amax',ytit='F770W',/xlog,yrange=[la_min(array_seds[*,3]),la_max(array_seds[*,3])]
  130 +
  131 +
  132 +window, 8,xsize=500,ysize=500
  133 +cgplot,amax,array_seds[*,5],psym='Filled Circle',color='black',xtit='amax',ytit='F1130W',/xlog,yrange=[la_min(array_seds[*,5]),la_max(array_seds[*,5])]
  134 +
  135 +ratio = array_seds[*,5]/array_seds[*,3]
  136 +
  137 +window, 9,xsize=500,ysize=500
  138 +cgplot,amax,ratio,psym='Filled Circle',color='black',xtit='amax',ytit='ratio',/xlog,yrange=[la_min(ratio),la_max(ratio)]
  139 +
  140 +
  141 +END
  142 +
  143 +
... ...