Commit 6768fb27fc73ea08bfd2cfe0765e64a6bc2334bc

Authored by Raphael Maris
1 parent 8a0c9c20
Exists in master

First comit

Showing 1 changed file with 173 additions and 0 deletions   Show diff stats
LabTools/IRAP/RM/phangs_make_herschel_images.pro 0 → 100644
... ... @@ -0,0 +1,173 @@
  1 +PRO phangs_make_herschel_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_Herschel_JWST_cycle1/'
  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 +name_loc_filters = ['PACS_BLUE','PACS_GREEN','PACS_RED','PSW']
  71 +herschel_filters=dustem_filter_names2filters(name_loc_filters)
  72 +Nherschel_filter=n_elements(herschel_filters)
  73 +herschel_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nherschel_filter)
  74 +
  75 +name_filters = ['pacs70','pacs100','pacs160','spire250']
  76 +
  77 +files=data_dir+use_source_name+'_scanamorphos_v25_'+name_filters+'_0.fits'
  78 +;NGC0628_F148_bkg_subtracted_mw_corrected.fits NGC0628_F169_bkg_subtracted_mw_corrected.fits NGC0628_N219_bkg_subtracted_mw_corrected.fits
  79 +;NGC0628_F154_bkg_subtracted_mw_corrected.fits NGC0628_F172_bkg_subtracted_mw_corrected.fits
  80 +wavs_mic=dustem_filter2wav(herschel_filters)
  81 +cmic=3.e14 ;mic/sec
  82 +
  83 +spire_beams_as = [450., 795., 1665.] ; in arcsec^2
  84 +
  85 +
  86 +omega_250 = spire_beams_as[0]*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
  87 +
  88 +omega_350 = spire_beams_as[1]*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
  89 +
  90 +omega_500 = spire_beams_as[2]*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
  91 +
  92 +imrange=[-0.5,2]*1e-17 ;erg/s/cm2/AA
  93 +imrange=[-0.5,1] ;MJy/sr
  94 +
  95 +FOR i=0L,Nherschel_filter-1 DO BEGIN
  96 + file=files[i]
  97 + d_narray=mrdfits(file,0,h)
  98 + d = d_narray[*,*,0]
  99 + d_err = d_narray[*,*,1]
  100 +
  101 +
  102 +
  103 +
  104 +
  105 +
  106 + ind=where(finite(d) NE 1,count)
  107 +
  108 + if (STRPOS('spire250', name_filters[i]) NE -1) then begin
  109 +
  110 + fact = 1.e-6/omega_250
  111 +
  112 +
  113 +
  114 + endif
  115 +
  116 + if (STRPOS('spire350', name_filters[i]) NE -1) then begin
  117 +
  118 + fact = 1.e-6/omega_350
  119 +
  120 + endif
  121 +
  122 + if (STRPOS('spire500', name_filters[i]) NE -1) then begin
  123 +
  124 + fact = 1/e-6/omega_500
  125 +
  126 + endif
  127 +
  128 + if (STRPOS(name_filters[i],'spire') EQ -1) then begin
  129 +
  130 + steradians=(sxpar(h,'CDELT2')/!radeg)^2 ;pixel in sr
  131 +
  132 + fact=1.e-6/steradians ;This is to rescale from /pixel to /sr
  133 +
  134 +
  135 +
  136 + endif
  137 +
  138 + IF count NE 0 THEN d[ind]=la_undef()
  139 + d=la_mul(d,fact)
  140 +
  141 + ;image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
  142 + IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
  143 + d=project2(h,d,href,/silent)
  144 + d_err=project2(h,d_err,href,/silent)
  145 + ENDIF
  146 + tit=use_source_name+' '+herschel_filters[i]
  147 + ;image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
  148 + herschel_images[*,*,0,i]=d
  149 + herschel_images[*,*,1,i]=d_err
  150 +
  151 +ENDFOR
  152 +
  153 +;Invent variances (will have to do better)
  154 +; perc_error=5./100.
  155 +; astrosat_images[*,*,1,*]=la_power(la_mul(astrosat_images[*,*,0,*],perc_error),2) ;assumed intensity variance
  156 +;=== check for null variances
  157 +ind=where(herschel_images[*,*,1,*] EQ 0.,count)
  158 +IF count NE 0 THEN BEGIN
  159 + message,'There are null variances ...',/continue
  160 + stop
  161 +ENDIF
  162 +
  163 +IF keyword_set(save) THEN BEGIN
  164 + file_save=save_data_dir+use_source_name+'_herschel_data.sav'
  165 + save,herschel_images,herschel_filters,href,file=file_save
  166 + message,'Saved '+file_save,/info
  167 +ENDIF
  168 +;
  169 +; mwrfits,herschel_images,save_data_dir+use_source_name+'_herschel_data.fits',href
  170 +
  171 +the_end:
  172 +
  173 +END
0 174 \ No newline at end of file
... ...