OV-GSO-DC / dustem-wrapper_idl

Blame view

LabTools/IRAP/JPB/fit_phangs_ngc0628_nir_continuum.pro 31.8 KB
Edit Raw Normal View History
d6343d56   Jean-Philippe Bernard   First commit 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
 
PRO fit_phangs_ngc0628_nir_continuum,from_restore=from_restore,brute_force=brute_force, $
									extract_seds=extract_seds,nostop=nostop,muse_stuff=muse_stuff,muse_data=muse_data,astrosat_data=astrosat_data, $
									show_plots=show_plots

;fit_phangs_ngc0628_nir_continuum,/nostop
;fit_phangs_ngc0628_nir_continuum,/muse_stuff
;fit_phangs_ngc0628_nir_continuum,/muse_data         ;do the Muse data "filter" save file
;fit_phangs_ngc0628_nir_continuum,/astrosat_data         ;do the Astrosat data save file
;fit_phangs_ngc0628_nir_continuum,/extract_seds
;fit_phangs_ngc0628_nir_continuum,/from_restore,/show_plots     ;fit seds of some voronoi bins
;fit_phangs_ngc0628_nir_continuum,/from_restore,/brute_force

;fit_phangs_ngc0628_nir_continuum,/muse_data

pdp_define_la_common
window,0
obp=[1.1,0,1.15,1]
source_name='ngc0628'
data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'

;stop

IF keyword_set(from_restore) THEN goto,from_restore
IF keyword_set(extract_seds) THEN goto,extract_seds
IF keyword_set(muse_stuff) THEN goto,muse_stuff
;IF keyword_set(extract_muse) THEN goto,NH_map
IF keyword_set(muse_data) THEN goto,muse_data
IF keyword_set(astrosat_data) THEN goto,astrosat_data

jwst_stuff:

;========== get JWST reference header
file=data_dir+source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits' 
d=mrdfits(file,1,h)
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
print,sxpar(h,'CDELT1')
href=cd2astro_header(h)
sxaddpar,href,'EQUINOX',2000.
image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
;stop
Nx=sxpar(href,'NAXIS1')
Ny=sxpar(href,'NAXIS2')

filters_names=['F200W','F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W']
filters=dustem_filter_names2filters(filters_names)
Nfilters=n_elements(filters)
jwst_images=fltarr(Nx,Ny,2,Nfilters)

IF not keyword_set(nostop) THEN stop

;========== JWST stuff
i=0L
file=data_dir+source_name+'_nircam_lv3_f200w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
help,d
;D               DOUBLE    = Array[7289, 9476]
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name

image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

;stop

file=data_dir+source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits' 
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
help,d
;D               DOUBLE    = Array[3515, 4576]
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_nircam_lv3_f335m_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_nircam_lv3_f360m_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f770w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.2,10],image_color_table='jpbloadct',/silent,tit=tit,off_bar=obp
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f1000w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f1130w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.3,20],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f2100w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

;WCO map:
file='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/ngc0628_12m+7m+tp_co21_broad_mom0.fits'
d=readfits(file,h)
;sxaddpar,h,'CTYPE1','RA---TAN'
;sxaddpar,h,'CTYPE2','RA---TAN'
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN WCO=project2(h,d,href,/silent) ELSE WCO=d
fact=4.e20/1.e21
NHCO=la_mul(WCO,fact)   ;NH from CO in 1e21 H/cm2
tit=source_name+' '+'NHCO [1e21 H/cm2]'
image_cont20,NHCO,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop


;Invent variances (will have to do better)
;jwst_images[*,*,1,*]=(jwst_images[*,*,0,*]*5./100.)^2   ;assumed intensity variance
perc_error=5./100.
jwst_images[*,*,1,*]=la_power(la_mul(jwst_images[*,*,0,*],perc_error),2)   ;assumed intensity variance
;=== check for null variances
ind=where(jwst_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
ENDIF

save,jwst_images,filters,href,NHCO,file=data_dir+'ngc0628_jwst_images.sav'

MUSE_data:

;=== this is just to get href
restore,data_dir+'ngc0628_jwst_images.sav',/verb

muse_data_dir='/Volumes/PILOT_FLIGHT1/PHANGS_MUSE/DR2p2/coopt/filterimages/NGC0628/'
;muse_filters=['COUSIN','DUPONT']
;muse_filters=['SDSS2','SDSS3','SDSS4']
muse_filters=dustem_filter_names2filters(['sdss_g','sdss_r','sdss_i'])
Nmuse_filter=n_elements(muse_filters)
muse_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nmuse_filter)
files=muse_data_dir+'NGC0628_'+['SDSS_gcopt','SDSS_rcopt','SDSS_icopt']+'.fits'
4c500bde   Jean-Philippe Bernard   updated with new ... 
223
224
225
226
 
;imrange=[-0.5,200]/2000.
imrange=[-1.e-4,0.001]
imrange=[-1.e-2,1.]
muse_undefined=0.
d6343d56   Jean-Philippe Bernard   First commit 
227
228
 
FOR i=0L,Nmuse_filter-1 DO BEGIN
	file=files[i]
4c500bde   Jean-Philippe Bernard   updated with new ... 
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
 
	;toto=mrdfits(file,0,hh)
	d=mrdfits(file,1,h)        ;here undefined values are at 0.
	hh=cd2astro_header(h)
	dv=mrdfits(file,2,hv)     ;is this looks like it's a variance
	mask=mrdfits(file,3,hm)
	ind=where(mask EQ 1,count)     ;this looks like a mask
	IF count NE 0 THEN BEGIN
		d[ind]=la_undef()
		dv[ind]=la_undef()
	ENDIF
	print,la_min(d),la_max(d)
	;     -7.37902      17186.9
	factor=dustem_muse_filter_conversion_factor(muse_filters[i],pivot_wav=pivot_wav)  ;goes from Muse filter units to muJy/pix
	omega_pix=(sxpar(hh,'CDELT2')/!radeg)^2    ;simeq 9e-13 sr
	;use_factor=factor/omega_pix*1.e-12         ;simeq 500
	;use_factor=factor/omega_pix*1.e-20         ;simeq 500
	use_factor=(1./factor)/omega_pix*1.e-12         ;simeq 500
	;jpb's estimate:
	cmic=2.99792458e14
	use_factor_jp=1.e-7*1.e4*pivot_wav*(pivot_wav/1.e4)/cmic/omega_pix
	print,factor,omega_pix,use_factor,use_factor_jp
	use_factor=use_factor_jp
	;stop
	;      492.255   9.4017732e-13   5.2357685e-06
    d=la_mul(d,use_factor)  ;Now supposidely in MJy/sr (but in what unit convention ??)
    dv=la_mul(dv,use_factor^2)  ;Now supposidely in (MJy/sr)^2
	print,la_min(d),la_max(d)
    ;-3.8634841e-05     0.089986754
    tit=source_name+' '+muse_filters[i]
	image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit,off_bar_pos=obp,axis_color_table=1
    ;stop
	IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
		d=project2(h,d,href,/silent,save_it='/tmp/save_project2.sav')
		dv=project2(h,dv,href,/silent,restore_it='/tmp/save_project2.sav')
	ENDIF
d6343d56   Jean-Philippe Bernard   First commit 
264
 
	muse_images[*,*,0,i]=d
4c500bde   Jean-Philippe Bernard   updated with new ... 
265
 
	muse_images[*,*,1,i]=dv
d6343d56   Jean-Philippe Bernard   First commit 
266
 
	;print,filter_name
4c500bde   Jean-Philippe Bernard   updated with new ... 
267
 
	image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit,off_bar_pos=obp,axis_color_table=1
d6343d56   Jean-Philippe Bernard   First commit 
268
269
270
271
 
ENDFOR

;Invent variances (will have to do better)
;CAUTION: some variances were 0 because some fluxes are 0, had to add abs_variance
4c500bde   Jean-Philippe Bernard   updated with new ... 
272
273
274
 
;perc_error=5./100.
;abs_variance=(1.e-12)^2
;muse_images[*,*,1,*]=la_add(la_power(la_mul(muse_images[*,*,0,*],perc_error),2),abs_variance)   ;assumed intensity variance
d6343d56   Jean-Philippe Bernard   First commit 
275
276
277
278
279
 
;=== check for null variances
ind=where(muse_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
4c500bde   Jean-Philippe Bernard   updated with new ... 
280
281
282
 
	toto=muse_images[*,*,1,*]
	toto[ind]=la_undef()
	muse_images[*,*,1,*]=toto
d6343d56   Jean-Philippe Bernard   First commit 
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
 
ENDIF

save,muse_images,muse_filters,href,file=data_dir+'ngc0628_muse_filters_data.sav'
message,'Saved '+data_dir+'ngc0628_muse_filters_data.sav',/info

stop

astrosat_data:
;BUNIT   = 'Angstrom-1 cm-2 erg s-1'
;in fact Angstrom-1 cm-2 erg s-1 PER PIXEL ?!

restore,data_dir+'ngc0628_jwst_images.sav',/verb

astrosat_data_dir='/Volumes/PILOT_FLIGHT1/PHANGS_ASTROSAT/v1p0/release/'
;muse_filters=['COUSIN','DUPONT']
astrosat_filters=dustem_filter_names2filters(['F148W','F154W','F169M','F172M','N219M'])
Nastrosat_filter=n_elements(astrosat_filters)
astrosat_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nastrosat_filter)
files=astrosat_data_dir+'NGC0628_'+['F148','F154','F169','F172','N219']+'_bkg_subtracted_mw_corrected.fits'
;NGC0628_F148_bkg_subtracted_mw_corrected.fits	NGC0628_F169_bkg_subtracted_mw_corrected.fits	NGC0628_N219_bkg_subtracted_mw_corrected.fits
;NGC0628_F154_bkg_subtracted_mw_corrected.fits	NGC0628_F172_bkg_subtracted_mw_corrected.fits
file_save_astrosat=data_dir+'ngc0628_astrosat_data.sav'
wavs_mic=dustem_filter2wav(astrosat_filters)
cmic=3.e14   ;mic/sec
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
imrange=[-0.5,2]*1e-17    ;erg/s/cm2/AA
imrange=[-0.5,1]   ;MJy/sr

FOR i=0L,Nastrosat_filter-1 DO BEGIN
	file=files[i]
	d=mrdfits(file,0,h)
	steradians=(sxpar(h,'CDELT2')/!radeg)^2   ;pixel in sr
	fact=facts[i]*4.*!pi/steradians
	ind=where(finite(d) NE 1,count)
	IF count NE 0 THEN d[ind]=la_undef()
	d=la_mul(d,fact)
	image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
	IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN d=project2(h,d,href,/silent)
	tit=source_name+' '+astrosat_filters[i]
	image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
	astrosat_images[*,*,0,i]=d
	stop
ENDFOR

;Invent variances (will have to do better)
perc_error=5./100.
astrosat_images[*,*,1,*]=la_power(la_mul(astrosat_images[*,*,0,*],perc_error),2)   ;assumed intensity variance
;=== check for null variances
ind=where(astrosat_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
ENDIF

save,astrosat_images,astrosat_filters,href,file=file_save_astrosat
message,'Saved '+file_save_astrosat,/info

stop


muse_stuff:

restore,data_dir+'ngc0628_jwst_images.sav',/verb

;========== MUSE stuff
;Muse stellar bins definition
;st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=NZ)
;print,Nage,NZ,Nbins
;          13           6          78
;st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,/young,Nbins=Nbins,Nage=Nage,NZ=NZ,bins=bins)
;print,Nage,NZ,Nbins
;          17           7         119
;print,la_min(bins),la_max(bins)
;            0          67
st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=NZ,bins=bins)
print,Nage,NZ,Nbins
;                    13           6          78
print,la_min(bins),la_max(bins)
;           0          77

;Muse weights
;st_muse_weights=read_muse_phangs_weights(object='NGC0628',/young,bin_numbers=bin_numbers)
;print,n_elements(bin_numbers)
;          68
st_muse_weights=read_muse_phangs_weights(object='NGC0628',bin_numbers=bin_numbers)
print,n_elements(bin_numbers)
;          78

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

image_cont20,voronoi_id,href,/square,/silent,image_color_table='jpbloadct',title='Voronoi Bins'

save,st_templates,st_muse_weights,voronoi_id,age_values,metalicity_values,bins,href,file=data_dir+'ngc0628_muse_images.sav'
message,'saved '+data_dir+'ngc0628_muse_images.sav',/info

stop

extract_seds:

restore,data_dir+'ngc0628_jwst_images.sav',/verb
restore,data_dir+'ngc0628_muse_images.sav',/verb
restore,data_dir+'ngc0628_astrosat_data.sav',/verb
restore,data_dir+'ngc0628_muse_filters_data.sav',/verb

stop

;=== extract and save observed seds in Muse Voronoi bins
extract_all_muse_phangs_seds,source_name,voronoi_id,jwst_images,filters,indices=all_seds_indices,file=data_dir+'ngc0628_jwst_seds_muse_pixels.sav',counts=counts
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb  ;just to get all_seds_indices
extract_all_muse_phangs_seds,source_name,voronoi_id,muse_images,muse_filters,use_these_indices=all_seds_indices,file=data_dir+'ngc0628_muse_seds_muse_pixels.sav',counts=counts
extract_all_muse_phangs_seds,source_name,voronoi_id,astrosat_images,astrosat_filters,use_these_indices=all_seds_indices,file=data_dir+'ngc0628_astrosat_seds_muse_pixels.sav',counts=counts

;==== The following is not really useful, can be done through SED aggregation later
;==== However, not equivalent to the above, because of undefined pixels in some of the images, and the way this is handled by sed extractor
stop
use_filters=[astrosat_filters,muse_filters,filters]
use_Nfilters=n_elements(use_filters)
use_images=fltarr([sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,use_Nfilters])
use_images[*,*,*,0:n_elements(filters)-1]=jwst_images
use_images[*,*,*,n_elements(filters):n_elements(filters)+n_elements(muse_filters)-1]=muse_images
use_images[*,*,*,n_elements(filters)+n_elements(muse_filters):n_elements(filters)+n_elements(muse_filters)+n_elements(astrosat_filters)-1]=astrosat_images

extract_all_muse_phangs_seds,source_name,voronoi_id,use_images,use_filters,indices=all_seds_indices,file=data_dir+'ngc0628_all_seds_muse_pixels.sav',counts=counts

stop

NH_map:
restore,data_dir+'ngc0628_jwst_images.sav',/verb
restore,data_dir+'ngc0628_muse_images.sav',/verb
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb

Nvor=max(voronoi_id)

ebv=st_muse_weights.reddening
Rv=3.1
Av_o_NH=1./2.  ;2 mag per 1e21 H/cm2
;Do an Av map
NH_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))  ;in 1e21 H/cm2
FOR vid=0LL,Nvor-1 DO BEGIN
	IF vid mod 100 EQ 0 THEN print,strtrim(vid/Nvor*100.,2)+' %'
	;ind=where(voronoi_id EQ vid)
    Nh_map[*all_seds_indices[vid]]=Rv*ebv[vid]/Av_o_NH
ENDFOR

save,NH_map,file=data_dir+'ngc0628_muse_NH.sav'

from_restore:

dustem_init,show_plots=show_plots
restore,data_dir+'ngc0628_jwst_images.sav',/verb
restore,data_dir+'ngc0628_muse_images.sav',/verb
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb
jwst_seds=all_seds
jwst_filters=(*jwst_seds[0]).filter
restore,data_dir+'ngc0628_muse_NH.sav',/verb
restore,data_dir+'ngc0628_astrosat_data.sav',/verb
restore,data_dir+'ngc0628_muse_seds_muse_pixels.sav',/verb
muse_seds=all_seds
muse_filters=(*muse_seds[0]).filter
restore,data_dir+'ngc0628_astrosat_seds_muse_pixels.sav',/verb
astrosat_seds=all_seds
astrosat_filters=(*astrosat_seds[0]).filter

Nvor=max(voronoi_id)
Nfilters=n_elements(jwst_filters)+n_elements(astrosat_filters)+n_elements(muse_filters)
all_seds=fltarr(Nfilters,Nvor)
4c500bde   Jean-Philippe Bernard   updated with new ... 
449
450
451
452
 
FOR i=0L,Nvor-1 DO BEGIN
	(*muse_seds[i]).stokesI=((*muse_seds[i]).stokesI)/242315.*1.e8
ENDFOR
FOR i=0L,Nvor-1 DO all_seds[*,i]=[(*astrosat_seds[i]).stokesI,((*muse_seds[i]).stokesI),(*jwst_seds[i]).stokesI]
d6343d56   Jean-Philippe Bernard   First commit 
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
 
;look for maximum astrosat value
ind=where(all_seds[2,*] EQ max(all_seds[2,*]))
print,ind
;       39665
print,all_seds[*,ind]
;      2.35202      2.75261      3.46727      3.68084      1.62485      1054.94      927.261      404.267      10.5454     -32768.0     -32768.0     -32768.0     -32768.0     -32768.0     -32768.0     -32768.0
ind=where(all_seds[2,*] GT 3.*la_mean(all_seds[2,*]))
print,ind[100]
;4699
print,all_seds[*,ind[100]]
;    0.0954950     0.126592     0.144323     0.120900     0.114967      68.6855      52.6373      38.7998      1.15606     0.638490     0.606459     0.529654      2.30200     0.782402      2.51761      1.86948
ind=where(all_seds[2,*]/all_seds[10,*] GT 1.)
print,ind[0]
;12306
print,all_seds[*,ind[0]]
;     0.323970     0.441198     0.591049     0.494629     0.182164      67.2941      50.9925      33.8923     0.967236     0.531406     0.498163     0.374723      1.39314     0.683776      1.90688      1.85371
print,ind[3]
print,all_seds[*,ind[3]]
;       12631
;     0.532320     0.654246     0.788243     0.643166     0.604815      85.0280      61.7483      41.8740      1.29023     0.754277     0.714327     0.569998      1.49024      1.07281      2.24375      1.84148
print,ind[4]
print,all_seds[*,ind[4]]
;       12640
;     0.421350     0.600934     0.639170     0.736419     0.392632      102.730      74.2435      50.7555      1.26155     0.704108     0.633040     0.518208      1.23460     0.844332      1.82023      2.09209
ind=where(all_seds[2,*]/all_seds[10,*] GT 0.6 and all_seds[2,*]/all_seds[10,*] LT 1.)
print,ind[0]
print,all_seds[*,ind[0]]
;        5960
;    0.0603282     0.103469     0.260063     0.171956    0.0811385      53.2903      43.1425      29.1814     0.689211     0.402946     0.368217     0.300235      1.35131     0.478370      1.73637      1.34104
print,ind[20]
print,all_seds[*,ind[20]]
;       13439
;     0.625768     0.544996     0.378714     0.838214     0.591592      92.5699      65.3618      41.9704      1.00230     0.566440     0.614708     0.448454      2.74491      1.15966      3.66431      2.25899


stop

IF keyword_set(brute_force) THEN goto,brute_force

;toto=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars,info_only=info_only,info_st=info_st)
;stop

;==== SFH ?? (Not sure what this is)
;file='/Volumes/PILOT_FLIGHT1/PHANGS_MUSE/DR2p2/coopt/MUSEDAP/fiducial/NGC0628-0.92asec_table_SFH.fits'
;toto=mrdfits(file,0,h) ;This is ?
;st_sfh=mrdfits(file,1,h) ;This is ?
;help,st_sfh,/str
;** Structure <478ef448>, 10 tags, length=72, data length=68, refs=1:
;   ID              LONG                 0
;   BIN_ID          LONG             47068
;   X               DOUBLE           121.80000
;   Y               DOUBLE          0.60000000
;   FLUX            DOUBLE           24.835750
;   SNR             DOUBLE           6.8757743
;   XBIN            DOUBLE           121.99930
;   YBIN            DOUBLE           2.9014085
;   SNRBIN          DOUBLE           103.86098
;   NSPAX           LONG               284
;print,minmax(st_sfh.id)
;           0        1425
;print,minmax(st_sfh.bin_id)
;       47028       47068


;vid=2000L ;selected voronoi ID
;vid=3000L ;selected voronoi ID
;vid=10000L
;vid=15000L ;selected voronoi ID
;vid=20000L ;selected voronoi ID
;vid=25000L
;vid=30000L ;selected voronoi ID
vid=35000L ;selected voronoi ID
use_NHmap=NHCO       ;used NH map in 1.e21 (from CO)

;galaxy_distance=10.  ;MPc

;use_model='DBP90'    ;Example with default keywords uses the DBP90 model
;use_model='MC10'    ;Example with default keywords uses the DBP90 model
use_model='DL07'    ;Example with default keywords uses the DBP90 model
use_polarization=0   ; initialize Dustemwrap in no polarization mode 
use_window=2          ; default graphics window number to use for plotting the results
use_verbose=0
use_Nitermax=5        ; maximum number of iterations for the fit
dustem_define_la_common
	;parameters to fit
pd = [ $
     '(*!dustem_params).G0', $                           ;G0
     '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$          ;PAH1 mass fraction
     ;'(*!dustem_params).grains(2).mdust_o_mh', $         ;VSGs
     ;'(*!dustem_params).grains(3).mdust_o_mh', $         ;VSGs
     ;'(*!dustem_params).grains(4).mdust_o_mh', $         ;BGs
     'dustem_plugin_phangs_stellar_continuum_1']         ;stellar continuum amplitude
;iv =   [1.6, [2.2e-3, 2.2e-3, 5.7e-3, 5.7e-3],2.e-3] 
iv =   [1.6, [2.2e-3, 2.2e-3],5.e-5] 
;iv =   [1.6, [2.2e-3, 2.2e-3,2e-3, 2e-3],5.e-5] 
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10]
llims=[1.e-2,[1.e-4,1.e-4],1.e-10]
;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10]

pd_to_update=['dustem_plugin_phangs_stellar_continuum_1','(*!dustem_params).grains(0).mdust_o_mh']
pd_filter_names=['NIRCAM11','MIRI2']

;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots
;stop
!dustem_nocatch=1
!dustem_verbose=use_verbose
IF keyword_set(noobj) THEN !dustem_noobj=1
!EXCEPT=2 ; for debugging
;=== INFORMATION TO MAKE THE PLOT
;yr=[1.00e-2,10.0] ; y-axis limits
yr=[1.00e-3,100.0] ; y-axis limits
;xr=[0.5,1.00e2] ; x-axis limits
xr=[0.1,1.00e2] ; x-axis limits
tit='FIT INTENSITY EXAMPLE' ; plot title
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
xtit=textoidl('\lambda (\mum)') ; x-axis title

;voronoi bins for test. Last series are randomly picked
vids=[5960,13439,12640,12631,12306,4699,[2000,3000,10000,15000,20000,25000,30000L,35000,400000,450000]]
Nvids=n_elements(vids)
one_sed=(*jwst_seds[0])[0]
Nfilters=n_elements(astrosat_filters)+n_elements(muse_filters)+n_elements(filters)
seds=replicate(one_sed,Nfilters)
FOR ii=0L,Nvids-1 DO BEGIN
	vid=vids[ii]
	;index=where(voronoi_id EQ vid,count)

	;===== get weights for the given voronoi bins
	weights=phangs_binid2weights(st_muse_weights,vid,st_templates,age_values,metalicity_values,reddening=reddening)
	;help,age_values,metalicity_values,weights,reddening

	;===== extract SED
	;There is no need to extract, since they already have bee extracted
	;sed=dustem_sed_extractor(jwst_images,index,filters,/total_intensity_only)
	Nf=0L & Nt=Nf+n_elements(astrosat_filters)-1
	seds[Nf:Nt]=*astrosat_seds[vid]
	Nf=Nt+1 & Nt=Nf+n_elements(muse_filters)-1
4c500bde   Jean-Philippe Bernard   updated with new ... 
595
596
 
	;seds[Nf:Nt]=*muse_seds[vid]
	seds[Nf:Nt]=(*muse_seds[vid])
d6343d56   Jean-Philippe Bernard   First commit 
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
 
	Nf=Nt+1 & Nt=Nf+n_elements(jwst_filters)-1
	seds[Nf:Nt]=*jwst_seds[vid]
	sed=seds
	;needed to add this because astrosat had 0 variances
	sed.sigmaii=(5./100.*sed.stokesI)^2
	;stop
	;===== fixed parameters
	fpd=['dustem_plugin_phangs_stellar_continuum_2']   ;This is Amplitude, E(B-V)
	fiv=[reddening] ;This is Aplitude, E(B-V)    ;WHY 2.e-4 ?
	Nage=n_elements(age_values)
	NZ=n_elements(metalicity_values)
	NH_value=la_mul(la_mean(use_NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2
	sed.stokesI=sed.stokesI/NH_value   ;normalize SED to 1.e20 H/cm2
	sed.sigmaII=sed.sigmaII/NH_value^2   ;normalize SED to 1.e20 H/cm2
	;stop
	FOR i=0L,Nage-1 DO BEGIN
		FOR j=0L,NZ-1 DO BEGIN
			IF weights[i,j] NE 0 AND bins[i,j] NE la_undef() THEN BEGIN
				ij=lonarr(1,2) & ij[0,0]=i & ij[0,1]=j
				ind=ij2index(ij,[Nage,Nz])+1
				fpd=[fpd,'dustem_plugin_phangs_stellar_continuum_'+strtrim(ind[0]+2,2)]   ;+2 is because amplitude,E(B-V) is parameter 0 of pluggin and bins start at zero
				fiv=[fiv,weights[i,j]]   
			ENDIF
		ENDFOR
	ENDFOR
    ;stop
	;== SET THE OBSERVATIONAL STRUCTURE
	dustem_set_data,m_fit=sed,m_show=sed
     ;stop
    ;get a first guess at parameters (must be done after dustem_set_data)
	new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol)
	stop
	;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
	;== ADJUSTED DURING THE FIT
	dustem_init_params,use_model,pd,new_iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims
    tit='Voronoi '+strtrim(vid,2)
	;===  RUN THE FIT
	;stop
	t1=systime(0,/sec)
	res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
                      ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
                      ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
                      ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plots)
	t2=systime(0,/sec)
    IF not keyword_set(nostop) THEN stop
ENDFOR

message,'========== Finished individual fitting',/continue
stop


;==== Brute force
brute_force:

Nvor=max(voronoi_id)

;ebv=st_muse_weights.reddening
;Rv=3.1
;Av_o_NH=1./2.  ;2 mag per 1e21 H/cm2
;;Do an Av map
;NH_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))  ;in 1e21 H/cm2
;FOR vid=0LL,Nvor-1 DO BEGIN
;	IF vid mod 100 EQ 0 THEN print,strtrim(vid/Nvor*100.,2)+' %'
;	ind=where(voronoi_id EQ vid)
;    Nh_map[ind]=Rv*ebv[vid]/Av_o_NH
;ENDFOR

;stop

dir=!dustem_wrap_soft_dir+'/Grids/'
table_name=dir+'DBP90_JWST_G0_YPAH_YVSG_4Phangs.fits'
;table_name=dir+'TEST_DBP90_JWST_G0_YPAH_YVSG_4Phangs.fits'
;Nvor=20000 ;for test

GOs=fltarr(Nvor)
Ypahs=fltarr(Nvor)
YVSGs=fltarr(Nvor) 
facts=fltarr(Nvor)
chi2s=fltarr(Nvor)
rchi2s=fltarr(Nvor)
use_NHmap=NHCO       ;used NH map in 1.e21 (from CO)
;use_NHmap=NH_map       ;used NH map in 1.e21 (from MUSE)

;filters_names=['F200W','F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W']
;filters=dustem_filter_names2filters(filters_names)
;==== filters used for brute force fit
use_filters_names=['F0770W','F1000W','F1130W','F2100W']
use_filters=dustem_filter_names2filters(use_filters_names)
ind_filters=[4,5,6,7]   ;to be done better

;show_sed=1
show_sed=0
FOR vid=0LL,Nvor-1 DO BEGIN
   IF vid mod 10 EQ 0 THEN BEGIN
   	 message,'Doing sed '+strtrim(vid,2)+' '+strtrim(1.*vid/Nvor*100.,2)+' %',/continue
   ENDIF
   sed=*all_seds[vid]
   sed=sed[ind_filters]
   ;=== normalize to NH=1.e20 H/cm2
   NH_value=la_mul(la_mean(use_NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2
   sed.stokesI=la_div(sed.stokesI,NH_value)
   sed.sigmaII=la_div(sed.sigmaII,NH_value^2)
   ;index=where(voronoi_id EQ vid,count)
   ;sed=dustem_sed_extractor(jwst_images,index,filters,/total_intensity_only)
   params=dustem_brute_force_fit(sed,table_name,use_filters,fact=fact,chi2=chi2,/normalize,rchi2=rchi2,show_sed=show_sed)
   GOs[vid]=params[0]
   Ypahs[vid]=params[1]/fact ;extenssive quantities must be devided by normalization factor
   Yvsgs[vid]=params[2]/fact
   facts[vid]=fact
   chi2s[vid]=chi2
   rchi2s[vid]=rchi2
   print,params[0],params[1],params[2],fact,chi2
ENDFOR

stop

;save,GOs,Ypahs,Yvsgs,file=dir+source_name+'DBP90_JWST_G0_YPAH_YVSG.fits'
save,GOs,Ypahs,Yvsgs,file=dir+source_name+'TEST_DBP90_JWST_G0_YPAH_YVSG.fits'

;stop
;=== make maps of parameters
chi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
rchi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
GO_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
Ypah_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
YVSG_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
fact_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
FOR vid=0LL,Nvor-1 DO BEGIN
	IF vid mod 100 EQ 0 THEN print,(1.*vid)/Nvor*100.
	;ind=where(voronoi_id EQ vid)
    GO_map[*all_seds_indices[vid]]=GOs[vid]
    Ypah_map[*all_seds_indices[vid]]=Ypahs[vid]
    YVSG_map[*all_seds_indices[vid]]=Yvsgs[vid]
    chi2_map[*all_seds_indices[vid]]=chi2s[vid]
    rchi2_map[*all_seds_indices[vid]]=rchi2s[vid]
    fact_map[*all_seds_indices[vid]]=facts[vid]
ENDFOR

;stop
jwst_4_coutours=fltarr((size(jwst_images))[1],(size(jwst_images))[2],2)
jwst_4_coutours[*,*,0]=jwst_images[*,*,0,5]
jwst_4_coutours[*,*,1]=jwst_images[*,*,0,5]
levs=create_contour_st(jwst_4_coutours,lev1=[1,5])
levs.(0).color=100
levs.(1).color=200
levs.(1).thick=2

obp=[1.1,0.,1.15,1]
image_cont20,jwst_4_coutours[*,*,0],href,/square,imrange=[-1.,10],image_color_table='jpbloadct',/silent,title='',off_bar=obp,axis_color_table=1,levels=levs

window,0,xsize=800,ysize=1000
aa=alog10(GO_map)
wset,0 & image_cont20,GO_map,href,/square,imrange=[-1.,100],image_color_table='jpbloadct',/silent,title='G0',off_bar=obp,axis_color_table=1
wset,0 & image_cont20,aa,href,/square,imrange=[0.5,2.5],image_color_table='jpbloadct',/silent,title='G0',off_bar=obp,axis_color_table=1,levels=levs
window,1,xsize=800,ysize=1000
wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.005,0.5],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1,levels=levs
window,2,xsize=800,ysize=1000
wset,2 & image_cont20,Yvsg_map,href,/square,imrange=[-0.005,0.5],image_color_table='jpbloadct',/silent,title='Yvsg',off_bar=obp,axis_color_table=1,levels=levs
window,3,xsize=800,ysize=1000
wset,3 & image_cont20,alog10(chi2_map),href,/square,imrange=[2,6],image_color_table='jpbloadct',/silent,title='chi2',off_bar=obp,axis_color_table=1,levels=levs
window,4,xsize=800,ysize=1000
image_cont20,alog10(fact_map),href,/square,imrange=[-1,5],image_color_table='jpbloadct',/silent,title='fact',off_bar=obp,axis_color_table=1,levels=levs

stop

wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.0001,0.005],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1
wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.0001,0.1],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1

wset,2 & image_cont20,Yvsg_map,href,/square,imrange=[-0.001,0.015],image_color_table='jpbloadct',/silent,title='Yvsg',off_bar=obp,axis_color_table=1

window,5
!p.multi=[0,3,2]
histGO=histogram(GOs,locations=GOv,Nbins=30)
cgplot,GOv,histGO,psym=10,/ylog,xtit='G0',yrange=[1,1.e5]
histYpah=histogram(Ypahs,locations=Ypahsv,Nbins=10)
cgplot,Ypahsv,histYpah,psym=10,/ylog,xtit='Ypah'
histYvsg=histogram(Yvsgs,locations=Yvsgsv,Nbins=10)
cgplot,Yvsgsv,histYvsg,psym=10,/ylog,xtit='Yvsg'
histchi2=histogram(chi2s,locations=chisv,Nbins=10)
cgplot,chisv,histchi2,psym=10,/ylog,xtit='Chi2'
histfacts=histogram(facts,locations=factsv,Nbins=10)
cgplot,factsv,histfacts,psym=10,/ylog,xtit='facts'


save,GO_map,Ypah_map,Yvsg_map,href,chi2_map,fact_map,GOs,Ypahs,Yvsgs,chi2s,facts,file=data_dir+'ngc0628_brute_force_fit_results.sav'

restore,data_dir+'ngc0628_brute_force_fit_results.sav',/verb


Nid=la_max(voronoi_id)
ids=fltarr(Nid)
snrs=fltarr(Nid)
counts=fltarr(Nid)
FOR i=0L,Nid-1 DO BEGIN
	ind=where(voronoi_id EQ i,count)
	ids[i]=i
	snrs[i]=snr[ind[0]]
	counts[i]=count
ENDFOR

cgplot,ids,snrs
cgplot,ids,counts
cgplot,counts,snrs

ind1=where(snr EQ max(snr))
ind=where(voronoi_id EQ voronoi_id[ind1],count)
print,count

ind=where(voronoi_id EQ 0,count)


file='/Volumes/PILOT_FLIGHT1/PHANGS_MUSE/DR2p2/coopt/MUSEDAP/fiducial/NGC0628-0.92asec_ppxf_SFH.fits'
st=mrdfits(file,1,h)
help,st,/str


END