DG Notes on 20170904 Event - Revision history

Dgary: /* Converting fits files to IDL save file */

2018-11-17T20:46:00Z

Converting fits files to IDL save file

← Older revision		Revision as of 20:46, 17 November 2018
Line 673:		Line 673:
	return		return
	end		end

			== Less Aggressive Restoring Beam ==
			I was not happy with the restoring beam I used when making the full suite of maps, which was the one Bin used for the 2017-Sep-10 event. Since I did like the images shown in the figures above, I have adjusted the parameters to be:
			sbeam = 40
			bm=max(sbeam*cfreqs[1]/cfreq,10.)
			In addition, I thought the clean box was a bit off, and maybe too large, so I changed that to:
			mask='circle [ [ 123pix ,118pix] , 45pix ]',

			I hope this will result in somewhat better images, and maybe fewer spurious sources in some maps.

Dgary: /* Converting fits files to IDL save file */

2018-11-12T15:55:29Z

Converting fits files to IDL save file

← Older revision		Revision as of 15:55, 12 November 2018
Line 620:		Line 620:
	; Get first map		; Get first map
	fits2map,files[0],omap		fits2map,files[0],omap
	~~sub_map~~,~~omap~~,rms~~,xrange=[400,500],yrange~~=[~~-160,-60~~]		rmsidx = where(shift(dist(256),128,128) gt 90)
	sdev = ~~stdev~~(rms~~.data~~)		rms = omap.data[rmsidx]
			sdev = sqrt((moment(rms,/nan))[1])
	sub_map,omap,map,xrange=[100,400],yrange=[-375,-75]		sub_map,omap,map,xrange=[100,400],yrange=[-375,-75]
	; Add various properties to this "template" map		; Add various properties to this "template" map
Line 655:		Line 656:
	eomaps[k-1,idx].id = 'EOVSA ' + string(freq,format='(F5.2)') + ' GHz'		eomaps[k-1,idx].id = 'EOVSA ' + string(freq,format='(F5.2)') + ' GHz'
	eomaps[k-1,idx].freq = freq		eomaps[k-1,idx].freq = freq
	~~sub_map,omap,~~rms~~,xrange=[400,500],yrange~~=[~~-160,-60~~]		rms = omap.data[rmsidx]
	sdev = ~~stdev~~(rms~~.data~~)		sdev = sqrt((moment(rms,/nan))[1])
	eomaps[k-1,idx].snr = max(map.data)/sdev		eomaps[k-1,idx].snr = max(map.data)/sdev
	; Get scale factor for this band and time		; Get scale factor for this band and time

Dgary: /* Converting fits files to IDL save file */

2018-11-12T13:11:21Z

Converting fits files to IDL save file

← Older revision		Revision as of 13:11, 12 November 2018
Line 635:		Line 635:
	add_prop,map,stokes='XX'		add_prop,map,stokes='XX'
	map.id = 'EOVSA 0.00 GHz'		map.id = 'EOVSA 0.00 GHz'
	; Allocate memory for all maps (30 frequencies and ~~150~~ times)		; Allocate memory for all maps (30 frequencies and 600 times)
	eomaps = replicate(map,30,600)		eomaps = replicate(map,30,600)
	k = 0		k = 0

Dgary: /* Converting fits files to IDL save file */

2018-11-12T13:10:26Z

Converting fits files to IDL save file

← Older revision		Revision as of 13:10, 12 November 2018
Line 605:		Line 605:
	# comparing the integrated flux density in the "used" part of the map with a total power spectrum ('eovsa_20170904_dynspec.sav') to obtain a scale factor to match the map flux with the total power spectrum,		# comparing the integrated flux density in the "used" part of the map with a total power spectrum ('eovsa_20170904_dynspec.sav') to obtain a scale factor to match the map flux with the total power spectrum,
	# writing the scale factor into each map structure		# writing the scale factor into each map structure
			Here is the complete procedure:
			pro fits2sav2step
			; Run in /srg/dgary/eovsa/solar/20170904_Mflare on baozi
			restore,'eovsa_20170904_dynspec.sav',/ver
			blist = fix((fghz - 2.86)*2)
			t_tp = anytim(etime)
			; Get file names in folder
			cd,'allimages2step'
			files = file_search('s01_*.fits',count=n)
			; Get time list of sorted filenames for band 1 (600 times with no gaps)
			tlist = strmid(files,13,6) ; format is 'HHMMSS'
			; Read all files (bands + times)
			files = file_search('*.fits',count=n)
			; Get first map
			fits2map,files[0],omap
			sub_map,omap,rms,xrange=[400,500],yrange=[-160,-60]
			sdev = stdev(rms.data)
			sub_map,omap,map,xrange=[100,400],yrange=[-375,-75]
			; Add various properties to this "template" map
			add_prop,map,rms=sdev
			add_prop,map,snr=max(map.data)/sdev
			add_prop,map,scalefac=1.0
			add_prop,map,rmsxran=[400.,500.]
			add_prop,map,rmsyran=[-160.,-60.]
			add_prop,map,dataunits='K'
			add_prop,map,rmsunits='K'
			add_prop,map,freq=2.92 + 0.5
			add_prop,map,frequnit='GHz'
			add_prop,map,stokes='XX'
			map.id = 'EOVSA 0.00 GHz'
			; Allocate memory for all maps (30 frequencies and 150 times)
			eomaps = replicate(map,30,600)
			k = 0
			; Loop over maps
			for i = 0, n-1 do begin
			reads,strmid(files[i],1,2),k ; Read band id as integer k
			fits2map,files[i],omap
			sub_map,omap,map,xrange=[100,400],yrange=[-375,-75]
			freq = 2.92 + k*0.5
			mapflux = 21.38e-1(freq1e9)^2/(3e8)^2 (2./206265.)^2 * total(map.data)
			idx = where(tlist eq strmid(files[i],13,6),cnt)
			if cnt eq 1 then begin
			print,k,idx
			eomaps[k-1,idx].data = map.data
			eomaps[k-1,idx].time = map.time
			eomaps[k-1,idx].l0 = map.l0
			eomaps[k-1,idx].b0 = map.b0
			eomaps[k-1,idx].rsun = map.rsun
			eomaps[k-1,idx].id = 'EOVSA ' + string(freq,format='(F5.2)') + ' GHz'
			eomaps[k-1,idx].freq = freq
			sub_map,omap,rms,xrange=[400,500],yrange=[-160,-60]
			sdev = stdev(rms.data)
			eomaps[k-1,idx].snr = max(map.data)/sdev
			; Get scale factor for this band and time
			tidx = (where(anytim(map.time) lt t_tp))[0]
			tpflux = mean(med[tidx,where(blist eq k)])
			fac = tpflux/mapflux
			; Write into map the scale factor needed to match TP
			eomaps[k-1,idx].data = map.data
			eomaps[k-1,idx].rms = sdev
			eomaps[k-1,idx].scalefac = fac
			endif
			endfor
			cd,'/srg/dgary/eovsa/solar/20170904_Mflare'
			save,file='20170904_2step.sav',eomaps
			return
			end

Dgary: /* Making the full set of images */

2018-11-12T13:09:42Z

Making the full set of images

← Older revision		Revision as of 13:09, 12 November 2018
Line 595:		Line 595:

	After the imaging is completed, the flux integrated over the map must be compared with the total power spectrum to find scaling factors to adjust the spectra. A single scale factor is determined for the entire flare for each frequency. This should result in smoother spectra for fitting.		After the imaging is completed, the flux integrated over the map must be compared with the total power spectrum to find scaling factors to adjust the spectra. A single scale factor is determined for the entire flare for each frequency. This should result in smoother spectra for fitting.

			== Converting fits files to IDL save file ==
			After all of the images are completed, I run the following IDL procedure to convert them to a single IDL save file. This IDL procedure features the following:
			# adding several keywords to the map structure,
			# calculating the RMS deviations in one corner of the map,
			# cutting out the central 300" x 300" portion of the map (151 x 151 pixels),
			# calculating the signal to noise ratio (snr) in the map, where the signal is defined as the maximum in the "used" part of the map,
			# putting the frequency, rms, and snr in each map structure,
			# comparing the integrated flux density in the "used" part of the map with a total power spectrum ('eovsa_20170904_dynspec.sav') to obtain a scale factor to match the map flux with the total power spectrum,
			# writing the scale factor into each map structure

Dgary: /* Second Stage */

2018-11-12T04:22:10Z

Second Stage

← Older revision		Revision as of 04:22, 12 November 2018
Line 509:		Line 509:

	The result after running this procedure is a quite pleasing improvement, as shown in these two summary image matrices.		The result after running this procedure is a quite pleasing improvement, as shown in these two summary image matrices.
	[[File:slf_t0_n0.png\|thumb\|~~600px~~\|Figure 1: Images in all 30 bands for one snapshot period (20:45:00~20:45:10 UT), after an iteration of phase-only selfcal.]]		[[File:slf_t0_n0.png\|thumb\|400px\|Figure 1: Images in all 30 bands for one snapshot period (20:45:00~20:45:10 UT), after an iteration of phase-only selfcal.]]
	[[File:slf_t0_n2.png\|thumb\|~~600px~~\|Figure 2: Same as Fig. 1, but after two additional amplitude selfcal iterations. The improvement is substantial, especially at the higher frequencies.]]		[[File:slf_t0_n2.png\|thumb\|400px\|Figure 2: Same as Fig. 1, but after two additional amplitude selfcal iterations. The improvement is substantial, especially at the higher frequencies.]]

	The final step is to apply the calibration tables to the original ms, to create the final selfcaled ms. The script section for this is		The final step is to apply the calibration tables to the original ms, to create the final selfcaled ms. The script section for this is

Dgary: /* Making the full set of images */

2018-11-12T04:20:56Z

Making the full set of images

← Older revision		Revision as of 04:20, 12 November 2018
Line 533:		Line 533:
	== Making the full set of images ==		== Making the full set of images ==

	The imaging step is straightforward, although time consuming since there are so many images to make (~~27000~~ images). The script is:		The imaging step is straightforward, although time consuming since there are so many images to make (18000 images). The script is:
	from suncasa.utils import helioimage2fits as hf		from suncasa.utils import helioimage2fits as hf
	import os		import os
	import numpy as np		import numpy as np
	from ptclean_cli import ptclean_cli as ptclean		from ptclean_cli import ptclean_cli as ptclean
	########### Note do applycal cannot be run with doptclean in the same sesssion ######!		########### Note do applycal cannot be run with
			############doptclean in the same sesssion ######!
	doptclean=1		doptclean=1
	workdir='/srg/dgary/eovsa/solar/20170904_Mflare/'		workdir='/srg/dgary/eovsa/solar/20170904_Mflare/'
Line 592:		Line 593:
	print 'cleaning spw '+sp+' unsuccessful. Proceed to next spw'		print 'cleaning spw '+sp+' unsuccessful. Proceed to next spw'
	continue		continue

			After the imaging is completed, the flux integrated over the map must be compared with the total power spectrum to find scaling factors to adjust the spectra. A single scale factor is determined for the entire flare for each frequency. This should result in smoother spectra for fitting.

Dgary: /* Making the full set of images */

2018-11-11T21:12:05Z

Making the full set of images

← Older revision		Revision as of 21:12, 11 November 2018
Line 532:		Line 532:

	== Making the full set of images ==		== Making the full set of images ==

			The imaging step is straightforward, although time consuming since there are so many images to make (27000 images). The script is:
			from suncasa.utils import helioimage2fits as hf
			import os
			import numpy as np
			from ptclean_cli import ptclean_cli as ptclean
			########### Note do applycal cannot be run with doptclean in the same sesssion ######!
			doptclean=1
			workdir='/srg/dgary/eovsa/solar/20170904_Mflare/'
			slfcaledms = workdir+'msdata/slfcal2/IDB20170904_2023_2103.ms.XX.peak.slfcaled.slfcaled2'
			spws=[str(s+1) for s in range(30)]
			antennas='0~10,12'
			trange=''
			pol='XX'
			if doptclean:
			imgdir=workdir+'/allimages2step/'
			spws=[str(s+1) for s in range(30)]
			tb.open(slfcaledms+'/SPECTRAL_WINDOW')
			reffreqs=tb.getcol('REF_FREQUENCY')
			bdwds=tb.getcol('TOTAL_BANDWIDTH')
			cfreqs=reffreqs+bdwds/2.
			tb.close()
			sbeam=30.
			for s,sp in enumerate(spws):
			cfreq=cfreqs[int(sp)]
			bm=max(sbeam*cfreqs[0]/cfreq,6.)
			imgprefix=imgdir+'s'+sp.zfill(2)+'_'
			print 'cleaning spw {0:s} with beam size {1:.1f}"'.format(sp,bm)
			try:
			ptclean(vis=slfcaledms,
			antenna=antennas,
			imageprefix=imgprefix,
			ncpu=20,
			twidth=4,
			doreg=True,
			usephacenter=False,
			toTb=True,
			spw=sp,
			mode='mfs',
			timerange=trange,
			imagermode='csclean',
			psfmode='clark',
			imsize=[256,256],
			cell=['2arcsec'],
			niter=1000,
			gain=0.05,
			stokes=pol,
			weighting='briggs',
			robust=0,
			restoringbeam=[str(bm)+'arcsec'],
			phasecenter='J2000 10h53m57 06d56m47',
			#mask='box [ [ 75pix , 90pix] , [205pix, 165pix ] ]',
			mask='circle [ [ 128pix ,128pix] , 50pix ]',
			pbcor=True,
			interactive=False,
			usescratch=False)
			except Exception as e:
			print e
			print 'cleaning spw '+sp+' unsuccessful. Proceed to next spw'
			continue

Dgary: /* Making the full set of images */

2018-11-11T20:42:45Z

Making the full set of images

← Older revision		Revision as of 20:42, 11 November 2018
Line 531:		Line 531:
	I made some images for other times from this selfcaled ms, and was at first a little disappointed, but I was using the clean masked created for the time 20:45:00 UT, which I now realize was too restricted in size. It nicely includes the double-footpoint source, but not the much weaker emission of the larger active region. The clean box needs to be enlarged to encompass this entire region, and in fact I can just go back to the large clean mask from stage 1, namely <code>mask='circle [[128pix ,128pix], 60pix]'</code>.		I made some images for other times from this selfcaled ms, and was at first a little disappointed, but I was using the clean masked created for the time 20:45:00 UT, which I now realize was too restricted in size. It nicely includes the double-footpoint source, but not the much weaker emission of the larger active region. The clean box needs to be enlarged to encompass this entire region, and in fact I can just go back to the large clean mask from stage 1, namely <code>mask='circle [[128pix ,128pix], 60pix]'</code>.

	=== Making the full set of images ===		== Making the full set of images ==

Dgary: /* Second Stage */

2018-11-11T20:32:30Z

Second Stage

← Older revision		Revision as of 20:32, 11 November 2018
Line 520:		Line 520:
	slftbs.append(tb[0])		slftbs.append(tb[0])
	clearcal(refms_slfcal)		clearcal(refms_slfcal)
	applycal(vis=refms_slfcal,gaintable=slftbs,spw=','.join(spws),selectdata=True,\		applycal(vis=refms_slfcal,gaintable=slftbs,
	antenna=antennas,interp='nearest',flagbackup=False,applymode='calonly',calwt=False)		spw=','.join(spws),selectdata=True,
			antenna=antennas,interp='nearest',
			flagbackup=False,applymode='calonly',
			calwt=False)
	if os.path.exists(refms_slfcaled):		if os.path.exists(refms_slfcaled):
	os.system('rm -rf '+refms_slfcaled)		os.system('rm -rf '+refms_slfcaled)