OV-GSO-DC / dustem-wrapper_idl

Commit 9beef1e99fd1aa9642f9dbdf49a7fec0afe70c5d

Authored by Jean-Philippe Bernard
1 parent 694026f6
Exists in master

first commit

Showing 3 changed files with 499 additions and 0 deletions   Show diff stats
Inline Side-by-side
LabTools/IRAP/JPB/phangs_make_astrosat_images.pro 0 → 100644
  Diff comments   View file @9beef1e
... ... @@ -0,0 +1,116 @@
  1 +PRO phangs_make_astrosat_images,source_name=source_name,save=save,show_images=show_images,nostop=nostop,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; phangs_make_astrosat_images
  6 +; PURPOSE:
  7 +; makes astrosat images usable by the Phangs ISRF project
  8 +; CATEGORY:
  9 +; Dustem Phangs
  10 +; CALLING SEQUENCE:
  11 +; phangs_make_astrosat_images[,source_name=][,/save][,/show_images][,/nostop]
  12 +; INPUTS:
  13 +; None
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; source_name : source name (default = 'ngc0628')
  16 +; OUTPUTS:
  17 +; None
  18 +; OPTIONAL OUTPUT PARAMETERS:
  19 +; None
  20 +; ACCEPTED KEY-WORDS:
  21 +; help = If set, print this help
  22 +; save = If set, save result
  23 +; show_images= if set, show images
  24 +; nostop = if set, does not stop
  25 +; COMMON BLOCKS:
  26 +; None
  27 +; SIDE EFFECTS:
  28 +; A file is written
  29 +; RESTRICTIONS:
  30 +; None
  31 +; PROCEDURE:
  32 +; input astrosat images are asumed to be in Angstrom-1 cm-2 erg s-1 PER PIXEL (even though BUNIT = 'Angstrom-1 cm-2 erg s-1' in header)
  33 +;
  34 +; EXAMPLES
  35 +; phangs_make_astrosat_images,source_name='ngc0628',/show,/nostop
  36 +; phangs_make_astrosat_images,source_name='ngc0628',/save,/show,/nostop
  37 +; MODIFICATION HISTORY:
  38 +; Written by J.-Ph. Bernard (2023)
  39 +; Evolution details on the DustEMWrap gitlab.
  40 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  41 +;-
  42 +
  43 +IF keyword_set(help) THEN BEGIN
  44 + doc_library,'phangs_make_astrosat_images'
  45 + goto,the_end
  46 +ENDIF
  47 +
  48 +use_source_name='ngc0628'
  49 +IF keyword_set(source_name) THEN use_source_name=source_name
  50 +
  51 +pdp_define_la_common
  52 +;window,0
  53 +obp=[1.1,0,1.15,1]
  54 +win=0L
  55 +
  56 +data_dir=!phangs_data_dir+'/phangs_drive/PHANGS_ASTROSAT/v1p0/release/'
  57 +save_data_dir=!phangs_data_dir+'/ISRF/WORK/'
  58 +
  59 +astrosat_data:
  60 +;BUNIT = 'Angstrom-1 cm-2 erg s-1'
  61 +;in fact Angstrom-1 cm-2 erg s-1 PER PIXEL ?!
  62 +
  63 +;This is to get the reference header
  64 +restore,save_data_dir+'ngc0628_jwst_images.sav',/verb
  65 +;% RESTORE: Restored variable: JWST_IMAGES.
  66 +;% RESTORE: Restored variable: FILTERS.
  67 +;% RESTORE: Restored variable: HREF.
  68 +;% RESTORE: Restored variable: NHCO.
  69 +
  70 +astrosat_filters=dustem_filter_names2filters(['F148W','F154W','F169M','F172M','N219M'])
  71 +Nastrosat_filter=n_elements(astrosat_filters)
  72 +astrosat_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nastrosat_filter)
  73 +files=data_dir+strupcase(use_source_name)+'_'+['F148','F154','F169','F172','N219']+'_bkg_subtracted_mw_corrected.fits'
  74 +;NGC0628_F148_bkg_subtracted_mw_corrected.fits NGC0628_F169_bkg_subtracted_mw_corrected.fits NGC0628_N219_bkg_subtracted_mw_corrected.fits
  75 +;NGC0628_F154_bkg_subtracted_mw_corrected.fits NGC0628_F172_bkg_subtracted_mw_corrected.fits
  76 +wavs_mic=dustem_filter2wav(astrosat_filters)
  77 +cmic=3.e14 ;mic/sec
  78 +facts=wavs_mic/cmic*(wavs_mic*1.e4)/4./!pi*1.e-7*1.e4*1.e20 ;from Flambda in ergs/cm2/s/AA to MJy/sr (note: was divided by 4*pi)
  79 +imrange=[-0.5,2]*1e-17 ;erg/s/cm2/AA
  80 +imrange=[-0.5,1] ;MJy/sr
  81 +
  82 +FOR i=0L,Nastrosat_filter-1 DO BEGIN
  83 + file=files[i]
  84 + d=mrdfits(file,0,h)
  85 + steradians=(sxpar(h,'CDELT2')/!radeg)^2 ;pixel in sr
  86 + fact=facts[i]*4.*!pi/steradians ;This is to rescale from /pixel to /sr
  87 + ind=where(finite(d) NE 1,count)
  88 + IF count NE 0 THEN d[ind]=la_undef()
  89 + d=la_mul(d,fact)
  90 + image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
  91 + IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN d=project2(h,d,href,/silent)
  92 + tit=source_name+' '+astrosat_filters[i]
  93 + image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
  94 + astrosat_images[*,*,0,i]=d
  95 + stop
  96 +ENDFOR
  97 +
  98 +;Invent variances (will have to do better)
  99 +perc_error=5./100.
  100 +astrosat_images[*,*,1,*]=la_power(la_mul(astrosat_images[*,*,0,*],perc_error),2) ;assumed intensity variance
  101 +;=== check for null variances
  102 +ind=where(astrosat_images[*,*,1,*] EQ 0.,count)
  103 +IF count NE 0 THEN BEGIN
  104 + message,'There are null variances ...',/continue
  105 + stop
  106 +ENDIF
  107 +
  108 +IF keyword_set(save) THEN BEGIN
  109 + file_save=save_data_dir+use_source_name+'_astrosat_data.sav'
  110 + save,astrosat_images,astrosat_filters,href,file=file_save
  111 + message,'Saved '+file_save,/info
  112 +ENDIF
  113 +
  114 +the_end:
  115 +
  116 +END
0 117 \ No newline at end of file
... ...
LabTools/IRAP/JPB/phangs_make_jwst_images.pro 0 → 100644
  Diff comments   View file @9beef1e
... ... @@ -0,0 +1,232 @@
  1 +PRO phangs_make_jwst_images,source_name=source_name,save=save,show_images=show_images,nostop=nostop
  2 +
  3 +;phangs_make_jwst_images,source_name='ngc0628',/save,/show,/nostop
  4 +
  5 +use_source_name='ngc0628'
  6 +IF keyword_set(source_name) THEN use_source_name=source_name
  7 +
  8 +pdp_define_la_common
  9 +;window,0
  10 +obp=[1.1,0,1.15,1]
  11 +win=0L
  12 +
  13 +;data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'
  14 +NH_data_dir=!phangs_data_dir+'/phangs_drive/PHANGS-JWST/DR1/'
  15 +data_dir=!phangs_data_dir+'/phangs_drive/PHANGS-JWST/DR1/'
  16 +save_data_dir=!phangs_data_dir+'/ISRF/WORK/'
  17 +
  18 +;========== get JWST reference header
  19 +file=data_dir+use_source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits'
  20 +d=mrdfits(file,1,h)
  21 +ind=where(finite(d) NE 1,count)
  22 +IF count NE 0 THEN d[ind]=la_undef()
  23 +print,sxpar(h,'CDELT1')
  24 +href=cd2astro_header(h)
  25 +sxaddpar,href,'EQUINOX',2000.
  26 +IF keyword_set(show_images) THEN BEGIN
  27 + window,win & win=win+1
  28 + image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
  29 +ENDIF
  30 +;stop
  31 +Nx=sxpar(href,'NAXIS1')
  32 +Ny=sxpar(href,'NAXIS2')
  33 +
  34 +filters_names=['F200W','F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W']
  35 +filters=dustem_filter_names2filters(filters_names)
  36 +Nfilters=n_elements(filters)
  37 +jwst_images=fltarr(Nx,Ny,2,Nfilters)
  38 +
  39 +IF not keyword_set(nostop) THEN stop
  40 +
  41 +;========== JWST stuff
  42 +i=0L
  43 +file=data_dir+use_source_name+'_nircam_lv3_f200w_i2d_anchor_atgauss1.fits'
  44 +d=mrdfits(file,0,h0)
  45 +d=mrdfits(file,1,h)
  46 +sxaddpar,h,'EQUINOX',2000.
  47 +ind=where(finite(d) NE 1,count)
  48 +IF count NE 0 THEN d[ind]=la_undef()
  49 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  50 +jwst_images[*,*,0,i]=d & i=i+1
  51 +help,d
  52 +;D DOUBLE = Array[7289, 9476]
  53 +print,sxpar(h,'EXTNAME')
  54 +filter_name=sxpar(h0,'FILTER')
  55 +print,filter_name
  56 +tit=source_name+' '+filter_name
  57 +IF keyword_set(show_images) THEN BEGIN
  58 + window,win & win=win+1
  59 + image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
  60 + IF not keyword_set(nostop) THEN stop
  61 +ENDIF
  62 +
  63 +;stop
  64 +
  65 +file=data_dir+use_source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits'
  66 +d=mrdfits(file,0,h0)
  67 +d=mrdfits(file,1,h)
  68 +sxaddpar,h,'EQUINOX',2000.
  69 +ind=where(finite(d) NE 1,count)
  70 +IF count NE 0 THEN d[ind]=la_undef()
  71 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  72 +jwst_images[*,*,0,i]=d & i=i+1
  73 +help,d
  74 +;D DOUBLE = Array[3515, 4576]
  75 +print,sxpar(h,'EXTNAME')
  76 +filter_name=sxpar(h0,'FILTER')
  77 +print,filter_name
  78 +tit=source_name+' '+filter_name
  79 +IF keyword_set(show_images) THEN BEGIN
  80 + window,win & win=win+1
  81 + image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
  82 + IF not keyword_set(nostop) THEN stop
  83 +ENDIF
  84 +
  85 +file=data_dir+use_source_name+'_nircam_lv3_f335m_i2d_anchor_atgauss1.fits'
  86 +d=mrdfits(file,0,h0)
  87 +d=mrdfits(file,1,h)
  88 +sxaddpar,h,'EQUINOX',2000.
  89 +ind=where(finite(d) NE 1,count)
  90 +IF count NE 0 THEN d[ind]=la_undef()
  91 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  92 +jwst_images[*,*,0,i]=d & i=i+1
  93 +print,sxpar(h,'EXTNAME')
  94 +filter_name=sxpar(h0,'FILTER')
  95 +print,filter_name
  96 +tit=source_name+' '+filter_name
  97 +IF keyword_set(show_images) THEN BEGIN
  98 + window,win & win=win+1
  99 + image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
  100 + IF not keyword_set(nostop) THEN stop
  101 +ENDIF
  102 +
  103 +file=data_dir+use_source_name+'_nircam_lv3_f360m_i2d_anchor_atgauss1.fits'
  104 +d=mrdfits(file,0,h0)
  105 +d=mrdfits(file,1,h)
  106 +sxaddpar,h,'EQUINOX',2000.
  107 +ind=where(finite(d) NE 1,count)
  108 +IF count NE 0 THEN d[ind]=la_undef()
  109 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  110 +jwst_images[*,*,0,i]=d & i=i+1
  111 +print,sxpar(h,'EXTNAME')
  112 +filter_name=sxpar(h0,'FILTER')
  113 +print,filter_name
  114 +tit=source_name+' '+filter_name
  115 +IF keyword_set(show_images) THEN BEGIN
  116 + window,win & win=win+1
  117 + image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
  118 + IF not keyword_set(nostop) THEN stop
  119 +ENDIF
  120 +
  121 +file=data_dir+use_source_name+'_miri_lv3_f770w_i2d_anchor_atgauss1.fits'
  122 +d=mrdfits(file,0,h0)
  123 +d=mrdfits(file,1,h)
  124 +sxaddpar,h,'EQUINOX',2000.
  125 +ind=where(finite(d) NE 1,count)
  126 +IF count NE 0 THEN d[ind]=la_undef()
  127 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  128 +jwst_images[*,*,0,i]=d & i=i+1
  129 +print,sxpar(h,'EXTNAME')
  130 +filter_name=sxpar(h0,'FILTER')
  131 +print,filter_name
  132 +tit=source_name+' '+filter_name
  133 +IF keyword_set(show_images) THEN BEGIN
  134 + window,win & win=win+1
  135 + image_cont20,d,href,/square,imrange=[-0.2,10],image_color_table='jpbloadct',/silent,tit=tit,off_bar=obp
  136 + IF not keyword_set(nostop) THEN stop
  137 +ENDIF
  138 +
  139 +file=data_dir+use_source_name+'_miri_lv3_f1000w_i2d_anchor_atgauss1.fits'
  140 +d=mrdfits(file,0,h0)
  141 +d=mrdfits(file,1,h)
  142 +sxaddpar,h,'EQUINOX',2000.
  143 +ind=where(finite(d) NE 1,count)
  144 +IF count NE 0 THEN d[ind]=la_undef()
  145 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  146 +jwst_images[*,*,0,i]=d & i=i+1
  147 +print,sxpar(h,'EXTNAME')
  148 +filter_name=sxpar(h0,'FILTER')
  149 +print,filter_name
  150 +tit=source_name+' '+filter_name
  151 +IF keyword_set(show_images) THEN BEGIN
  152 + window,win & win=win+1
  153 + image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
  154 + IF not keyword_set(nostop) THEN stop
  155 +ENDIF
  156 +
  157 +file=data_dir+use_source_name+'_miri_lv3_f1130w_i2d_anchor_atgauss1.fits'
  158 +d=mrdfits(file,0,h0)
  159 +d=mrdfits(file,1,h)
  160 +sxaddpar,h,'EQUINOX',2000.
  161 +ind=where(finite(d) NE 1,count)
  162 +IF count NE 0 THEN d[ind]=la_undef()
  163 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  164 +jwst_images[*,*,0,i]=d & i=i+1
  165 +print,sxpar(h,'EXTNAME')
  166 +filter_name=sxpar(h0,'FILTER')
  167 +print,filter_name
  168 +tit=source_name+' '+filter_name
  169 +IF keyword_set(show_images) THEN BEGIN
  170 + window,win & win=win+1
  171 + image_cont20,d,href,/square,imrange=[-0.3,20],image_color_table='jpbloadct',/silent,tit=tit
  172 + IF not keyword_set(nostop) THEN stop
  173 +ENDIF
  174 +IF not keyword_set(nostop) THEN stop
  175 +
  176 +file=data_dir+use_source_name+'_miri_lv3_f2100w_i2d_anchor_atgauss1.fits'
  177 +d=mrdfits(file,0,h0)
  178 +d=mrdfits(file,1,h)
  179 +sxaddpar,h,'EQUINOX',2000.
  180 +ind=where(finite(d) NE 1,count)
  181 +IF count NE 0 THEN d[ind]=la_undef()
  182 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
  183 +jwst_images[*,*,0,i]=d & i=i+1
  184 +print,sxpar(h,'EXTNAME')
  185 +filter_name=sxpar(h0,'FILTER')
  186 +print,filter_name
  187 +tit=source_name+' '+filter_name
  188 +IF keyword_set(show_images) THEN BEGIN
  189 + window,win & win=win+1
  190 + image_cont20,d,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
  191 + IF not keyword_set(nostop) THEN stop
  192 +ENDIF
  193 +IF not keyword_set(nostop) THEN stop
  194 +
  195 +;Invent variances (will have to do better)
  196 +;jwst_images[*,*,1,*]=(jwst_images[*,*,0,*]*5./100.)^2 ;assumed intensity variance
  197 +perc_error=5./100.
  198 +jwst_images[*,*,1,*]=la_power(la_mul(jwst_images[*,*,0,*],perc_error),2) ;assumed intensity variance
  199 +;=== check for null variances
  200 +ind=where(jwst_images[*,*,1,*] EQ 0.,count)
  201 +IF count NE 0 THEN BEGIN
  202 + message,'There are null variances ...',/continue
  203 + stop
  204 +ENDIF
  205 +
  206 +;WCO map:
  207 +file=NH_data_dir+use_source_name+'_12m+7m+tp_co21_broad_mom0.fits'
  208 +;file='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/'+use_source_name+'_12m+7m+tp_co21_broad_mom0.fits'
  209 +d=readfits(file,h)
  210 +;sxaddpar,h,'CTYPE1','RA---TAN'
  211 +;sxaddpar,h,'CTYPE2','RA---TAN'
  212 +sxaddpar,h,'EQUINOX',2000.
  213 +ind=where(finite(d) NE 1,count)
  214 +IF count NE 0 THEN d[ind]=la_undef()
  215 +IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN WCO=project2(h,d,href,/silent) ELSE WCO=d
  216 +fact=4.e20/1.e21
  217 +NHCO=la_mul(WCO,fact) ;NH from CO in 1e21 H/cm2
  218 +tit=source_name+' '+'NHCO [1e21 H/cm2]'
  219 +IF keyword_set(show_images) THEN BEGIN
  220 + window,win & win=win+1
  221 + image_cont20,NHCO,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
  222 + IF not keyword_set(nostop) THEN stop
  223 +ENDIF
  224 +
  225 +IF keyword_set(save) THEN BEGIN
  226 + save_file=save_data_dir+use_source_name+'_jwst_images.sav'
  227 + save,jwst_images,filters,href,NHCO,file=save_file
  228 + message,'Saved '+save_file,/continue
  229 +ENDIF
  230 +
  231 +
  232 +END
0 233 \ No newline at end of file
... ...
LabTools/IRAP/JPB/phangs_make_muse_filters_images.pro 0 → 100644
  Diff comments   View file @9beef1e
... ... @@ -0,0 +1,151 @@
  1 +PRO phangs_make_muse_filters_images,source_name=source_name,save=save,show_images=show_images,nostop=nostop
  2 +
  3 +;phangs_make_muse_filters_images,source_name='ngc0628',/show,/nostop
  4 +;phangs_make_muse_filters_images,source_name='ngc0628',/save,/show,/nostop
  5 +
  6 +use_source_name='ngc0628'
  7 +IF keyword_set(source_name) THEN use_source_name=source_name
  8 +
  9 +pdp_define_la_common
  10 +;window,0
  11 +obp=[1.1,0,1.15,1]
  12 +win=0L
  13 +
  14 +;data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'
  15 +data_dir=!phangs_data_dir+'/phangs_drive/PHANGS_MUSE/DR2p2/coopt/filterimages/NGC0628/'
  16 +save_data_dir=!phangs_data_dir+'/ISRF/WORK/'
  17 +
  18 +MUSE_data:
  19 +
  20 +;=== this is just to get href
  21 +restore,save_data_dir+'ngc0628_jwst_images.sav',/verb
  22 +
  23 +;muse_data_dir='/Volumes/PILOT_FLIGHT1/PHANGS_MUSE/DR2p2/coopt/filterimages/NGC0628/'
  24 +;muse_filters=['COUSIN','DUPONT']
  25 +;muse_filters=['SDSS2','SDSS3','SDSS4']
  26 +muse_filters=dustem_filter_names2filters(['sdss_g','sdss_r','sdss_i'])
  27 +Nmuse_filter=n_elements(muse_filters)
  28 +muse_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nmuse_filter)
  29 +files=data_dir+'NGC0628_'+['SDSS_gcopt','SDSS_rcopt','SDSS_icopt']+'.fits'
  30 +
  31 +;imrange=[-0.5,200]/2000.
  32 +imrange=[-1.e-4,0.001]
  33 +imrange=[-1.e-2,1.]
  34 +muse_undefined=0.
  35 +FOR i=0L,Nmuse_filter-1 DO BEGIN
  36 + file=files[i]
  37 + ;toto=mrdfits(file,0,hh)
  38 + d=mrdfits(file,1,h) ;here undefined values are at 0.
  39 + hh=cd2astro_header(h)
  40 + dv=mrdfits(file,2,hv) ;is this looks like it's a variance
  41 + mask=mrdfits(file,3,hm)
  42 + ind=where(mask EQ 1,count) ;this looks like a mask
  43 + IF count NE 0 THEN BEGIN
  44 + d[ind]=la_undef()
  45 + dv[ind]=la_undef()
  46 + ENDIF
  47 + print,la_min(d),la_max(d)
  48 + ; -7.37902 17186.9
  49 + factor=dustem_muse_filter_conversion_factor(muse_filters[i],pivot_wav=pivot_wav) ;goes from Muse filter units to muJy/pix
  50 + omega_pix=(sxpar(hh,'CDELT2')/!radeg)^2 ;simeq 9e-13 sr
  51 + ;use_factor=factor/omega_pix*1.e-12 ;simeq 500
  52 + ;use_factor=factor/omega_pix*1.e-20 ;simeq 500
  53 + use_factor=(1./factor)/omega_pix*1.e-12 ;simeq 500
  54 + ;jpb's estimate:
  55 + cmic=2.99792458e14
  56 + use_factor_jp=1.e-7*1.e4*pivot_wav*(pivot_wav/1.e4)/cmic/omega_pix
  57 + print,factor,omega_pix,use_factor,use_factor_jp
  58 + use_factor=use_factor_jp
  59 + ;stop
  60 + ; 492.255 9.4017732e-13 5.2357685e-06
  61 + d=la_mul(d,use_factor) ;Now supposidely in MJy/sr (but in what unit convention ??)
  62 + dv=la_mul(dv,use_factor^2) ;Now supposidely in (MJy/sr)^2
  63 + print,la_min(d),la_max(d)
  64 + ;-3.8634841e-05 0.089986754
  65 + tit=source_name+' '+muse_filters[i]
  66 + image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit,off_bar_pos=obp,axis_color_table=1
  67 + ;stop
  68 + IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
  69 + d=project2(h,d,href,/silent,save_it='/tmp/save_project2.sav')
  70 + dv=project2(h,dv,href,/silent,restore_it='/tmp/save_project2.sav')
  71 + ENDIF
  72 + muse_images[*,*,0,i]=d
  73 + muse_images[*,*,1,i]=dv
  74 + ;print,filter_name
  75 + IF keyword_set(show_images) THEN BEGIN
  76 + window,win & win=win+1
  77 + image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit,off_bar_pos=obp,axis_color_table=1
  78 + ENDIF
  79 +ENDFOR
  80 +
  81 +;Invent variances (will have to do better)
  82 +;CAUTION: some variances were 0 because some fluxes are 0, had to add abs_variance
  83 +;perc_error=5./100.
  84 +;abs_variance=(1.e-12)^2
  85 +;muse_images[*,*,1,*]=la_add(la_power(la_mul(muse_images[*,*,0,*],perc_error),2),abs_variance) ;assumed intensity variance
  86 +;=== check for null variances
  87 +ind=where(muse_images[*,*,1,*] EQ 0.,count)
  88 +IF count NE 0 THEN BEGIN
  89 + message,'There are null variances ...',/continue
  90 + stop
  91 + toto=muse_images[*,*,1,*]
  92 + toto[ind]=la_undef()
  93 + muse_images[*,*,1,*]=toto
  94 +ENDIF
  95 +
  96 +IF keyword_set(save) THEN BEGIN
  97 + file_save=save_data_dir+use_source_name+'_muse_filters_data.sav'
  98 + save,muse_images,muse_filters,href,file=file_save
  99 + message,'Saved '+file_save,/info
  100 +ENDIF
  101 +
  102 +muse_stuff:
  103 +
  104 +restore,save_data_dir+use_source_name+'_jwst_images.sav',/verb
  105 +
  106 +;========== MUSE stuff
  107 +st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=NZ,bins=bins)
  108 +print,Nage,NZ,Nbins
  109 +; 13 6 78
  110 +print,la_min(bins),la_max(bins)
  111 +; 0 77
  112 +;Muse weights
  113 +st_muse_weights=read_muse_phangs_weights(object='NGC0628',bin_numbers=bin_numbers)
  114 +print,n_elements(bin_numbers)
  115 +; 78
  116 +
  117 +voronoi_id=read_muse_phangs_voronoi_bins('NGC0628',header_in=href,snr_bin=snr_bin,snr_flux=snr_flux,flux=muse_flux) ;This is the voronoi bin map
  118 +
  119 +IF keyword_set(show_images) THEN BEGIN
  120 + window,win & win=win+1
  121 + image_cont20,voronoi_id,href,/square,/silent,image_color_table='jpbloadct',title='Voronoi Bins'
  122 +ENDIF
  123 +
  124 +IF keyword_set(save) THEN BEGIN
  125 + file_save=save_data_dir+use_source_name+'_muse_images.sav'
  126 + save,st_templates,st_muse_weights,voronoi_id,age_values,metalicity_values,bins,href,file=file_save
  127 + message,'saved '+file_save,/info
  128 +ENDIF
  129 +
  130 +NH_map:
  131 +
  132 +restore,save_data_dir+use_source_name+'_jwst_images.sav',/verb
  133 +restore,save_data_dir+use_source_name+'_muse_images.sav',/verb
  134 +restore,save_data_dir+use_source_name+'_jwst_seds_muse_pixels.sav',/verb
  135 +
  136 +Nvor=max(voronoi_id)
  137 +
  138 +ebv=st_muse_weights.reddening
  139 +Rv=3.1
  140 +Av_o_NH=1./2. ;2 mag per 1e21 H/cm2
  141 +;Do an Av map
  142 +NH_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2')) ;in 1e21 H/cm2
  143 +FOR vid=0LL,Nvor-1 DO BEGIN
  144 + IF vid mod 100 EQ 0 THEN print,strtrim(vid/Nvor*100.,2)+' %'
  145 + ;ind=where(voronoi_id EQ vid)
  146 + Nh_map[*all_seds_indices[vid]]=Rv*ebv[vid]/Av_o_NH
  147 +ENDFOR
  148 +
  149 +save,NH_map,file=save_data_dir+use_source_name+'_muse_NH.sav'
  150 +
  151 +END
0 152 \ No newline at end of file
... ...