PRO make_phangs_isrf_classes,bidon,source_name=source_name,fits=fits,help=help ;+ ; NAME: ; make_phangs_isrf_classes ; CALLING SEQUENCE: ; make_phangs_isrf_classes[,/save] ; PURPOSE: ; generates the ISRF classes used in fiting PHANGs data ; INPUTS: ; None ; OPTIONAL KEYWORDS: ; save = if set, save the classes into a .sav file ; fits = if set, saves the classes into a .fits file ; help = if set, print this help ; OUTPUTS: ; None ; OPTIONAL INPUT: ; None ; OPTIONAL OUTPUT: ; None ; PROCEDURE AND SUBROUTINE USED ; None ; SIDE EFFECTS: ; None ; EXAMPLE: ; make_phangs_isrf_classes,/save ; MODIFICATION HISTORY: ; written by Jean-Philippe Bernard ;- IF keyword_set(help) THEN BEGIN doc_library,'make_phangs_isrf_classes' goto,the_end ENDIF win=0L ;data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/' data_dir=!phangs_data_dir+'/ISRF/WORK/' use_model='DBP90' ;Example with default keywords uses the DBP90 model use_polarization=0 ; initialize Dustemwrap in no polarization mode ;== INITIALISE DUSTEM dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots use_source_name='ngc0628' IF keyword_set(source_name) THEN use_source_name=source_name ;==== restore needed data restore,data_dir+use_source_name+'_isrf_min_prediction.sav',/verb ;% RESTORE: Restored variable: ISRFS. ;% RESTORE: Restored variable: OBJECT_DISTANCE. ;% RESTORE: Restored variable: OBJECT_THICKNESS. ;% RESTORE: Restored variable: SOURCE_NAME. ;=== The following is restored only for variables HREF and ALL_SEDS_INDICES ;=== in order to generate the fits class map fo that galaxy restore,data_dir+use_source_name+'_ref_header.sav',/verb ;% RESTORE: Restored variable: HREF. ;=== The following is restored only for variables HREF and ALL_SEDS_INDICES ;=== in order to generate the fits class map fo that galaxy restore,data_dir+use_source_name+'_astrosat_seds_muse_pixels.sav',/verb ;% RESTORE: Restored variable: ALL_SEDS. ;% RESTORE: Restored variable: ALL_SEDS_INDICES. class_map=lonarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2')) ;stop ;=== Defines ISRF classes, based on a given galaxy lambir=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths) Nlamb=(size(ISRFS))[1] Nvor=(size(ISRFS))[2] ;Define normalized ISRFs n_ISRFS=fltarr(Nlamb,Nvor) ;ISRFs normalized to 1 micron n_wavs=fltarr(Nlamb,Nvor) ind=where(isrf_wavelengths GT 1.) indd=ind[0] FOR vid=0L,Nvor-1 DO BEGIN n_ISRFS[*,vid]=ISRFS[*,vid]/ISRFS[indd,vid] n_wavs[*,vid]=isrf_wavelengths ENDFOR ;=== Plot 2D histogram of normalized ISRFs in NGC0628 x=reform(n_wavs,Nlamb*Nvor) y=reform(n_ISRFS,Nlamb*Nvor) ;jpbloadct loadct,13 Nx=200 & Ny=200 Nx=100 & Ny=100 xmin=-1.5 & xmax=2. ymin=-5 & ymax=1. xtit='log10(wav [mic])' ytit='log10(ISRF/ISRF(1 mic))' range=[-1,5] window,win & win=win+1 plot_bin_correl,alog10(x),alog10(y),Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,xtit=xtit,ytit=ytit,range=range;,/xlog ;,/ylog ;=== compute the median ISRFS ;med_isrf=la_median(n_ISRFS,dim=-1) ;refs_wavs=[0.11] ;This is set to be just before the Lymann break reference_wavelength=1. ;all ISRFs are normalized to 1 mic refs_wavs=[0.15] ;This is set to be just before the Lymann break ;ref_ratios=interpol(med_isrf,isrf_wavelengths,refs_wavs) ;actually not used ;These are the ratios which will be used to classify ISRFs ;==== compute ISRF ration for each Voronoi bin isrf_ratio_values=fltarr(Nvor) FOR i=0L,Nvor-1 DO BEGIN isrf_ratio_values[i]=interpol(n_ISRFS[*,i],isrf_wavelengths,refs_wavs[0]) ENDFOR undef_array=fltarr(Nlamb) undef_array[*]=la_undef() one_class={number:0L $ ;class number ,name:'' $ ;class name ,ISRF_wavelengths:undef_array $ ;isrf wavelengths ,ISRF_template:undef_array $ ;isrf template, normalized at reference_wavelength ,flux_ratio:0. $ ;flux ratio between flux_ratio_wav and reference_wavelength ,flux_ratio_wav:0. $ ;flux ratio wavelength ,reference_wavelength:reference_wavelength $ ;ISRF normalisation wavelength ,flux_ratio_min:0. $ ;minimum flux ratio to belong to this class ,flux_ratio_max:0. $ ;maximum flux ratio to belong to this class } ; ,flux_ref:0. $ ; ; ,Nvor:0L $ ;number of Voronoi bins in the class ; ,voronoi_members:ptr_new() $ ;voronoi pixel members of that class ; ,voronoi_chi2:ptr_new()} ;chi2 between ref ISRF and voronoi member's ISRF class_ratio_range=[1.e-5,10.] ;reange of flux ratios used in the classification Nclasses=30L ;Number of classes used for classification ratios=range_gen(Nclasses,class_ratio_range,/log,/get_minmax,min_values=ratios_min,max_values=ratios_max) classes=replicate(one_class,Nclasses) ;====== fill in class values Nvor_class=lonarr(Nclasses) classes.number=lindgen(Nclasses) classes.flux_ratio=ratios classes.flux_ratio_wav=refs_wavs[0] classes.flux_ratio_min=ratios_min classes.flux_ratio_max=ratios_max median_isrf=fltarr(Nlamb) FOR i=0L,Nclasses-1 DO BEGIN ind=where(isrf_ratio_values GE ratios_min[i] AND isrf_ratio_values LT ratios_max[i],count) IF count NE 0 THEN BEGIN FOR j=0L,Nlamb-1 DO BEGIN median_isrf[j]=la_median(n_ISRFs[j,ind]) ENDFOR classes[i].ISRF_template=median_isrf classes[i].ISRF_wavelengths=isrf_wavelengths Nvor_class[i]=count FOR j=0L,count-1 DO BEGIN class_map[*ALL_SEDS_INDICES[ind[j]]]=i ENDFOR ENDIF ENDFOR ;stop ;=== Do a class (class 0) with a Mathis field class0=one_class class0.number=0L class0.name='Mathis' class0.reference_wavelength=1. class0.flux_ratio_wav=refs_wavs[0] file=!dustem_soft_dir+'/data/ISRF_MATHIS.DAT' ;file='/Users/jpb/Soft_Libraries/dustem_fortran/data/ISRF.DAT' readcol,file,Mathis_wavs,Mathis_ISRF Mathis_ISRF=interpol(Mathis_ISRF,Mathis_wavs,isrf_wavelengths) n_Mathis_ISRF=Mathis_ISRF/interpol(Mathis_ISRF,isrf_wavelengths,class0.reference_wavelength) ;class0.ISRF_template=ptr_new(n_Mathis_ISRF) class0.ISRF_template=n_Mathis_ISRF class0.ISRF_wavelengths=isrf_wavelengths ratio=interpol(n_Mathis_ISRF,isrf_wavelengths,class0.flux_ratio_wav) class0.flux_ratio=ratio class0.flux_ratio_min=ratio*(1.-0.1) ;arbitrary +-10% around ratio class0.flux_ratio_max=ratio*(1.+0.1) ;arbitrary +-10% around ratio Mathis_1mic=interpol(Mathis_ISRF,isrf_wavelengths,1.0) message,'Mathis (G0=1) has 4*pi*Inu='+strtrim(Mathis_1mic,2)+' ergs/s/cm2/Hz',/continue message,'Mathis (G0=1) has Inu(0.15 mic)/Inu(1 mic)='+strtrim(ratio,2)+' ',/continue ;=== Do a plot of ISRF classes window,win & win=win+1 colors=range_gen(Nclasses,[0,255]) yrange=[1.e-5,10] xrange=[1.e-2,1.e4] xtit='Wavlenghts [mic]' ytit='Normalized ISRF (Inu)' cgplot,class0.ISRF_wavelengths,class0.ISRF_template,/xlog,/ylog,yrange=yrange,xrange=xrange,xtit=xtit,ytit=ytit FOR i=0L,Nclasses-1 DO BEGIN oplot,classes[i].ISRF_wavelengths,classes[i].ISRF_template,color=colors[i] ENDFOR ;=== plot an histogram of the classes: window,win & win=win+1 xtit='ISRF class #' ytit='Number of Voronoi bins' cgplot,Nvor_class,psym=10,xtit=xtit,ytit=ytit ;=== plot the class map: cleanplot window,win & win=win+1 obp=[1.1,0.,1.15,1] imrange=[0,30] image_cont20,class_map,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar=obp,imrange=imrange ;stop ;save empty classes defined above ;IF keyword_set(save) THEN BEGIN ; file_save=data_dir+'isrf_classes_one_ratio.sav' ; file_save=data_dir+'isrf_classes_one_ratio_on_'+use_galaxy+'.sav' ; message,'Will save '+file_save,/continue ; message,'.c to proceed ',/continue ;stop ; save,classes,class0,file=file_save,/verb ; message,'Saved '+file_save,/continue ;ENDIF ;=== save a fits file IF keyword_set(fits) THEN BEGIN file_save=data_dir+'isrf_classes_one_ratio_on_'+use_source_name+'.fits' ;message,'Will save '+file_save,/continue ;message,'.c to proceed ',/continue ;stop mwrfits,class0,file_save,/create mwrfits,classes,file_save message,'Saved '+file_save,/continue file_save=data_dir+use_source_name+'_isrf_classes_map.fits' writefits,file_save,class_map,href message,'Saved '+file_save,/continue ENDIF the_end: END